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ABSTRACT 

We present numerical simulations of axisymmetric, magnetically driven outflows that 
reproduce the inferred properties of ultra-relativistic gamma-ray burst (GRB) jets. 
These results extend our previous simulations (Komissarov et al. 2007) of outflows ac- 
celerated to moderately relativistic speeds, which we applied to jets of active galactic 
nuclei. In contrast to several recent investigations, which have employed the magne- 
todynamics approximation, our numerical scheme solves the full set of equations of 
special-relativistic, ideal MHD, which enables us to explicitly calculate the jet velocity 
and magnetic-to-kinetic energy conversion efliciency — key parameters of interest for 
astrophysical applications. We conflrm that the magnetic acceleration scheme remains 
robust into the ultra-relativistic regime, as previously indicated by semi- analytic self- 
similar solutions. We find that all current-carrying outflows exhibit self-coUimation 
and consequent acceleration near the rotation axis, but that unconfined outflows lose 
causal connectivity across the jet and therefore do not coUimate or accelerate effi- 
ciently in their outer regions. We show that magnetically accelerated jets confined by 
an external pressure that varies as z~" (0 < a < 2) assume a paraboloidal shape 
z oc (where r, z are cylindrical coordinates and a > 1), and we obtain analytic 
expressions for the one-to-one correspondence between the pressure distribution and 
the asymptotic jet shape. We demonstrate that the acceleration efliciency of jets with 
paraboloidal streamlines is ^50%, with the numerical value being higher the lower 
the initial magnetization. We derive asymptotic analytic expressions for the acceler- 
ation of initially cold outflows along paraboloidal streamlines and verify that they 
provide good descriptions of the simulated flows. Our modelled jets (corresponding to 
3/2 < a < 3) attain Lorentz factors r>10^ on scales ~ 10^° — 10^^ cm, consistent 
with the possibility that long/soft GRB jets are accelerated within envelopes of col- 
lapsing massive stars, and F > 30 on scales ~ 9 x 10^ — 3 x 10^° cm, consistent with the 

possibility that short/hard GRB jets are accelerated on scales where they can be con- 
fined by moderately relativistic winds from accretion discs. We also find that F^v 1 
for magnetically accelerated jets, where 9^ is the half-opening angle of the poloidal 
streamlines, which implies that the 7-ray emitting components of GRB outflows are 
very narrow, with 9^ < 1° in regions where F > 100, and that the afterglow light curves 
of these components would either exhibit a very early jet break or show no jet break 
at all. 
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In the "standard" model of long-duration, soft-spectrum 
gamma-ray bursts (GRBs; e.g. |Piran| [2005[ ), the prompt 
high-energy emission arises in ultra-relativistic (bulk 
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Lorentz factor r>10^), highly colhmated (opening half- 
angle of a few degrees) jets. The high Lorentz factors are 
inferred from the requirement of a sufficiently low opacity 
to photon-photon annihilation or to scattering by photon 



annihilation-produced electron-positron pairs (e.g. Lithwick 
fc Sari|2001 1, whereas the jet opening angle is deduced from 
the detection of a panchromatic break in the light curve of 
the lower-energy afterglow emission (e.g. Rhoads 1999, Sari 
et al.|1999 |. Recent observations by the Swift satellite have 
indicated that various aspects of this model may need to be 
modified (e.g. Meszarosl|2006 Panaitescu|2007 Liang et al 
2008), but the basic picture of a colhmated F > 10 outflow 
IS still the accepted paradigm. 

Observations of long/soft GRBs and their afterglows 
have revealed that these events typically involve the release 
of a few times lO^'^ erg, although the fraction of this energy 
that corresponds to the 7-ray emitting outflow component 
may vary from source to source (e.g. [Berger et al.|[2003a[ 



Frail et al. 20051. The outflows in these GRBs have been 



argued to originate either in a magnetar or in a rapidly ac- 
creting stellar-mass black hole, formed in the collapse of a 
massive star. The jets could tap into the rotational energy 
of the neutron star, black hole or accretion disc through 
the agency of an ordered magnetic fleld that threads the 

^ 1994 Mcszaros & Rees 
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source (e.g. Usov| 1992 Thompson 
Katz" 1997 



Kluzniak & Rudorman 1998, Vlahaki s fc 
'Konigl.2001, 2003a b: Blandford.2002, Drenkhahn fc Spruit 



2002 Vlahakis at al. p003| Proga et al. p()03| McKinney 



2006 



2007 



Lyutikov|2006b Levinson'2006'; 'Ko missarov fc Barkov 
Bucciantini et al. 2008 ; Barkov fc Komissarov|2008 1 



For typical burst energies and durations the field ampli- 



tudes should be 



10^ 



10 G. Early models have pos- 



tulated that GRB outflows are driven purely thermally via 
annihilation of neutrinos emitted by the accretion disc. Al- 
though this model remains very popular, some recent studies 
have indicated that the neutrino heating may not be as ef- 



flcient as previously thought (e.g. Di Matteo et al. 2002 1 



At present, both the magnetic and the thermal mechanisms 
seem equally possible and it may well be that in many cases 
they operate simultaneously. In particular, neutrino heat- 
ing may play an important role in the initial acceleration of 
magnetized outflows (e.g. Vlahakis fc Konigl|2003a I and in 
determining their mass load (e.g. Levinson|2006 Barzilay fc 
Levinson|2008l ). 

While short/hard GRBs evidently have different pro- 
genitors (quite possibly merging neutron stars or neutron 
star/black hole pairs) and on average involve a smaller en- 
ergy release, a lower Lorentz factor, and weaker colhmation 
than long/soft GRBs, they may well represent the same ba- 
sic phenomenon and arise in relativistic outflows that are 
driven in a similar way (e.g. |Nakar|2007l ). 

The magnetic acceleration and colhmation of GRB out- 
flows needs to be studied within the framework of rela- 
tivistic magnetohydrodynamics (MHD). Although general- 
relativistic effects may influence the conditions near the base 
of the flow, most of the action takes place sufficiently far 
away from the central mass that the simpler equations of 
special-relativistic MHD can be employed. Since our focus 
in this paper is on the global structure of GRB jets, we 
henceforth consider only the special-relativistic theory. How- 
ever, even in this case there are qualitatively new effects in 
comparison with Newtonian MHD. These include the fact 



that, when the bulk Lorentz factor becomes large, the elec- 
tric force can no longer be neglected relative to the magnetic 
force and, in fact, becomes comparable to it in magnitude. 
Correspondingly, one needs to retain the displacement cur- 
rent and the electric charge density in Maxwell's equations. 
Another consequence of relativistic motion (which also af- 
fects unmagnetized flows) is the coupling between different 
spatial components of the momentum conservation equa- 
tion brought about by the appearance of the Lorentz fac- 
tor (which is calculated from the total velocity) in each of 
the component equations. Furthermore, in cases where the 
temperature (i.e. the characteristic velocity of internal mo- 
tions) is relativistic, one needs to take into account the en- 
thalpy contribution to the inertia of the flow. On account 
of these various factors, relativistic MHD does not naturally 
yield simple generalizations of results obtained in Newto- 
nian MHD. To simplify the treatment, various authors have 
adopted the force-free electrodynamics (also termed "mag- 
netodynamics" ) approximation, in which the matter inertia 
is neglected altogether. While this approach has led to useful 
insights and interesting exact solutions, it is inherently lim- 
ited in that one cannot explicitly calculate the fluid velocity 
and hence the efflciency of transforming electromagnetic en- 
ergy into kinetic form, which are key parameters of interest 
for astrophysical applications. 



In a pioneering work, Li et al. ( 1992 see also Con- 



topoulos||1994 1 derived exact semi-analjrtic MHD solutions 



of steady, axisymmetric, "cold" relativistic flows patterned 
after the Newtonian radially self-similar outflow solutions 



of Blandford fc Payne (19821. In contrast with the Newto- 



nian solutions, one cannot match the flow in the relativistic 
case to a given power-law radial distribution of the rotation 
velocity of the source (e.g. the oc r~^^^ rotation law of a 
Keplerian accretion disc) because the relativistic equations 
already contain the (constant) speed of light c. However, this 
constraint only affects the base of the flow (and, as shown 
by [Vlahakis fc Konigl|2003a[ it is possible to approximate a 
Keplerian disc even in this case by judiciously parametrizing 
the disc height above the origin of the coordinate system), 
and one can proceed to obtain the global structure of the 
outflow as in the Newtonian case. Li et al. ( 1992 1 identi- 



fled as a key property of the relativistic outflow solutions 
the spatially extended nature of the acceleration region, 
which continues well beyond the classical fast-magnetosonic 
surface. These results were further generalized to initially 
"hot" outflows by [Vlahakis fc Konigl| ( ["2003a|b[ ), who went 
on to apply the relativistic self-similar solutions to GRB out- 



flows (see also Vlahakis fc Konigl 2001 and Vlahakis at al 
|2003[ ) and to the lower-F jets imaged in active galactic nuclei 



(AGNs; [Vlahakis fc K6nigl|200"4t . The solutions obtained in 
these papers conflrmed that spatially extended acceleration 
is a generic property of MHD outflows that distinguishes 
it from purely hydrodynamic, thermally driven winds. [Vla-[ 



hakis fc Konigl (2001 2003a I noted that this property can 



be understood from the fact that the magnetic accelera- 
tion is determined from the joint solution of the Bernoulli 
equation (derived from the momentum conservation equa- 
tion along the poloidal magnetic fleld) and the trans-field 
equation (which describes the force balance in the trans- 
verse direction). The effective singular surface (the "event 
horizon" for the propagation of fast-magnetosonic waves) 
is the so-called modified fast magnetosonic surface, which 



Magnetic acceleration of ultra-relativistic jets in gamma-ray hurst sources 3 



can lie well beyond the corresponding classical surface. (The 
classical fast-magnetosonic surface is singular only when one 
solves the Bernoulli equation alone, assuming that the shape 
of the field lines is given; in Section [5] we further elaborate 
on the strong connection between acceleration and poloidal 
field-line shape in magnetically driven flows.) 

The semi-analytic solutions have also established the 
collimation properties of MHD outflows, demonstrating that 
they converge asymptotically to cylinders for flows that are 
Poynting flux-dominated at the source and to cones when 



the enthalpy flux is initially dominant (e.g. Vlahakis & 



[Konig l 2003a b). These solutions are, however, limited by 
the self-similarity assumption, which, besides restricting the 
angular velocity distribution at their base, also requires the 
magnetic flux distribution to be a power law in radius and 
only enables one current regime (current-carrying or return- 
current, but not a global current circuit) to be modelled 
by any given solution. To validate the applicability of these 
results under more realistic circumstances and to ascertain 
their dynamical stability, one needs to resort to numerical 
simulations. However, the large spatial extent of the accel- 
eration region (which, according to the semi-analytic solu- 
tions, typically covers several decades in spherical radius) 
has posed a strong challenge for such calculations: in fact, 
early attempts to simulate such flows were limited by numer- 
ical dissipation to maximum Lorentz factors that were only 
a small fraction (less than 1%) of the potentially achievable 
terminal value. 

[Komissarov et ah] ( |2007| hereafter Paper I) have taken 
a major step toward overcoming this challenge by employing 
a special-relativistic, ideal-MHD numerical scheme that was 
specifically designed to optimize accuracy and resolution and 
to minimize numerical dissipation. A key element of their ap- 
proach was the implementation of a grid-extension method 
that made it possible to follow the flow up to six decades 
in spatial scale while reducing the computation time by up 
to three orders of magnitude. They were able to model cold 
flows that converted nearly 80% of the initial Poynting flux 
into kinetic energy of Too ^ 10 baryons and demonstrated 
that the results were consistent with the available data on 
the acceleration of relativistic jets in AGNs. They found 
that the numerical solutions assumed a quasi-static configu- 
ration that was qualitatively in accord with the self-similar 
AGN jet models oflVlahakis & Konigll (l2004|. The simu- 



lations were, however, able to examine various aspects of 
the flow that could not be studied within the framework of 
a self-similar model (including the structure of outflows in 
which both the current and the return current flow within 
the jet and the dependence of the collimation properties on 
the shape of the jet boundary) and uncovered new features 
(such as the formation of a cylindrical core around the jet 
axis) that were inherently non-self-similar. 

In this paper we further extend the scheme presented 
in Paper I to cover the regime of GRB outflows. In particu- 
lar, we present simulations of outflows that attain terminal 
Lorentz factors Vac ^ 10^, following them over up to eight 
decades in axial scale. Besides cold jets, we also consider 
the case of an outfiow in which the enthalpy flux is a sig- 
nificant fraction of the injected energy flux. Owing to the 
larger range in F in comparison with the solutions presented 
in Paper I, the magnetic acceleration region can now be bet- 
ter isolated, which enables us to more accurately compare 



its behaviour with that of the self-similar solutions and to 
analyse it using the asymptotic forms of the Bernoulli and 
trans-field equations. We begin by reviewing the relativistic 
MHD formalism (Section[2| and our numerical scheme (Sec- 
tion [sj. We present key simulation results in Section |4] and 
discuss them in the context of the theory of magnetic accel- 
eration in Section[5] Section [6] deals with applications of our 
results to GRBs. Our conclusions are given in Section [7| 



2 BASIC EQUATIONS 

Since most of the acceleration takes place far away from the 
source, we assume that the space-time is flat. In an inertial 
frame at rest relative to the source, the relativistic ideal- 
MHD equations that describe the flow take the following 
form; continuity equation 



-gpu 



0, 



(1) 



where p is the rest mass density of matter, it"^ is its 4- 
velocity, and g is the determinant of the metric tensor; 
energy-momentum equations 



(l/c)a(^n) + d,{V^gT\) = 



-d„{gaB)T' 



(2) 



where T""' is the total stress-energy-momentum tensor; in- 
duction equation 



(l/c)a(B')-f e"^'*a,(i50 = 0, 



(3) 



where Cijk = ^/^eijk is the Levi-Civita tensor of the absolute 
space (ei23 = 1 for right-handed systems and £123 = — 1 for 
left-handed ones) and 7 is the determinant of the spatial part 
of the metric tensor {ytj = gij); the solenoidal condition 



0. 



(4) 



The total stress-energy-momentum tensor, T'^'^ , is a 
sum of the stress-energy momentum tensor of matter. 



(m) =wu u jc +pg 



(5) 



where p is the thermodynamic pressure and w is the enthalpy 
per unit volume, and the stress-energy momentum tensor of 
the electromagnetic field. 



1 

47r 



(F"^F„/5)/ 



(6) 



where F"'^ is the Maxwell tensor. The electric and magnetic 
fields are defined as measured by an observer stationary rel- 
ative to the spatial grid, which gives 



B' = 



-e^^'^Fi. 



and 

E^ = Fi . 

In the limit of ideal MHD 

Ei = -CijkV^ /c, 



(7) 
(8) 
(9) 



where = /u^ is the usual 3- velocity of the plasma. 

In all of our simulations we use an isentropic equation 
of state 



Qp" , 



(10) 
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where Q =const and s = 4/3. This relation enables us to ex- 
clude the energy equation from the integrated system. How- 
ever, the momentum equation remains intact, including the 
non-linear advection term. Therefore, if the conditions for 
shock formation were to arise, our calculation would cap- 
ture that shockj^The enthalpy per unit volume is 



2 , 
■ pC + 



(11) 



is the Poynting flux per unit rest-mass energy flux, and cr 
is the ratio of the Poynting flux to the hydrodynamic (rest- 
mass plus thermal plus kinetic) energy flux. For cold flows 
Q = 0, w — pc? . (Here and in the rest of the paper we use a 
hat symbol over vector indices to indicate their components 



in a normalized coordinate basis.) From equation (f6l it 
follows that the Lorentz factor F cannot exceed p. 



2.1 Field-line constants 

The poloidal magnetic field is fully described by the az- 
imuthal component of the vector potential, 



1 



ij't' 



(12) 



For axisymmetric solutions A^p — 'I'/27r, where *I'(a;'), the 
so-called magnetic flux function, is the total magnetic flux 
enclosed by the circle a;' =const (a;' being the coordinates 
of the meridional plane). Stationary and axisymmetric ideal 
MHD flows have five quantities that propagate unchanged 
along the magnetic field lines and thus are functions of 'I' 
alone. These are k, the mass flux per unit magnetic flux; 
fl, the angular velocity of magnetic field lines; /, the total 
angular momentum flux per unit rest-mass flux; p., the to- 
tal energy flux per unit rest-mass energy fiux; and Q, the 
entropy per particle: 



Bp 

r r B- 



I = 



p 
w 

+ r- 
2nkc pc^ 

p = Ph + Pni, 

and 

Q = P/P, 



(13) 

(14) 

(15) 
(16) 

(17) 



where Up = Tvp is the magnitude of the poloidal compo- 
nent of the 4-velocity, Bp is the magnitude of the poloidal 
component of the magnetic field, r is the cylindrical radius, 



(18) 



is the total electric current flowing through a loop of radius 
r around the rotation axis. 



w 
pc^ 



(19) 



is the total hydrodynamic energy (rest mass plus thermal 
plus kinetic) flux per unit rest-mass energy flux. 



Pn 



■ Phcy 



2-Kkc? 



(20) 



^ Since entropy is fixed, the compression of our shocks would be 
the same as for continuous compression waves. This would give 
a higher jump in density for the same jump in pressure than 
in a proper (dissipative) shock. Fortunately, we do not need to 
contend with this issue in practice as shocks do not form in our 
simulations. 



3 NUMERICAL SIMULATIONS 

To maintain a firm control over the jet's confinement and to 
prevent complications related to numerical diffusion of the 
denser plasma from the jet's surroundings, we study outflows 
that propagate inside a solid funnel of a prescribed shape[^ 
Specifically, we consider axisymmetric funnels 



where z and r are the cylindrical coordinates of the funnel 
wall and a — 2/3, 1, 3/2, 2 and 3. We employ elliptical 
coordinates {^,rj,(j}}, where 



and 



V 



-l/a 



r 2 

— + 2 

a 



(21) 



(22) 



(see Paper 1 for details). 

We use a Godunov-type numerical code based on the 



scheme described in Komissarov ( 1999a I. To reduce numeri- 



cal diffusion we apply parabolic reconstruction instead of the 
linear one of the original code. Our procedure, in brief, is to 
calculate minmod-averaged flrst and second derivatives and 
use the first three terms of the Taylor expansion for spatial 
reconstruction. This simple procedure results in a notice- 
able improvement in the solution accuracy even though the 
new scheme is still not 3rd-order accurate because of the 
non-uniformity of the grid. 

The grid is uniform in the £, direction (the polar angle 
direction when we use spherical coordinates) , where in most 
runs it has a total of 60 cells. To check the convergence, some 
runs were repeated with a doubled resolution. The cells are 
elongated in the rj direction (the radial direction when we 
use spherical coordinates), reflecting the elongation of the 
funnel. Very elongated cells lead to a numerical instability, 
so we imposed an upper limit of 40 on the length/width 
ratio. 

To speed up the simulations, we implement a section- 
ing of the computational grid as described in |Komissarov| 
& Lyubarsky (2004). In each section, which is shaped as a 



ring, the numerical solution is advanced using a time step 
based on the local Courant condition. It is twice as large as 
the time step of the adjacent inner ring and twice as small 
as the time step of the adjacent outer ring. This approach is 
particularly effective for conical flows but less so for highly 
collimated, almost cylindrical configurations. 

^ As was already noted in Paper I, in real astrophysical systems 
the shape of the boundary is determined by the spatial distri- 
bution of the pressure or the density of the confining ambient 
medium. The effective ambient pressure distributions implied by 
the adopted funnel shapes are considered in Section 15.31 



Magnetic acceleration of ultra-relativistic jets in gamma-ray hurst sources 5 



The equations are dimensionalized in the following man- 
ner. The unit of length, L, is such that 77^ = 1, where the 
subscript i refers to the inlet boundary. The unit of time is 
T = L/c. The unit of mass is M = L^Bo/47rc^, where Bo is 
the dimensional magnitude of the r] component of magnetic 
field at the inlet (so the dimensionless magnitude of at 
the inlet is \/47r). In applications, L is the typical length- 
scale of the launch region, T is the light crossing time of 
that region and Bo is the typical strength of the poloidal 
magnetic field at the origin. Notice that L does not have to 
be the size of the rotating object at the base of the jet and in 
particular it cannot be identified with the radius of the black 
hole event horizon which allows only inflows. When dimen- 
sional estimates are required we use the expected magnitude 
of the light-cylinder radius, ric = c/Q,. The mass scale M 
does not represent the mass of the central object but rather 
the rest-mass equivalent of the magnetic energy within the 
magnetosphere. 



3.1 Boundary conditions 

3.1.1 Inlet boundary 

We treat the inlet boundary, 7;^ = 1, as a surface of a per- 
fectly conducting rotator with either a uniform angular ve- 
locity = f^o or with 



n = f7o(l + a2(e/C,)'+a3(C/e.)'), 



(23) 



where the subscript j refers to the jet boundary (funnel 
wall). In this paper we set a2 = 0.778 and 03 = —1.778. 
The angular velocity profile is directly related to the distri- 
bution of the electric current in the jet, which for r ^ ric is 
given by 

1 



-f7B„r 



(24) 



(see Paper I, or equation 34 in Section 5.21. In fact, the 
current is driven by the electric fleld associated with the ro- 
tating poloidal field, and the electric charge conservation 
requires the circuit to eventually close. In the case of a 
constant the return current flows over the jet boundary, 
whereas in the case of differential rotation with f2(^j) = it 
flows mainly inside the jet (within 0.75 < ^/^j < 1 for the 



f2 distribution given by equation 231. The solid-body rota- 



tion law provides a very good description of the behaviour of 
magnetic field lines that thread the horizon of a black hole 
or the surface of a magnetized star. This choice is therefore 
entirely appropriate for the black-hole or magnetar theory 
of GRB jets. On the other hand, differential rotation is a 
natural choice for jets that are launched from an accretion 
disc, and although the distribution ( |23[ ) does not correspond 
to a realistic disc model, it should nevertheless capture the 
qualitative aspects of such a system]^ 

The condition of perfect conductivity allows us to fix the 
azimuthal component of the electric field and the r) compo- 
nent of the magnetic field: 



^ Note in tills connection tliat [Tchekhovskoy et al.| ( |2008| l simu- 
lated a force- free black- liole/disc outflow in which current flowed 
out along field lines that threaded the uniformly rotating hole 
and returned along field lines attached to the differentially rotat- 
ing disc. 



= 0, B'' = Bo at 77 = 77, . 

From the first of these conditions we derive 

,fi 



V 



V 



and (using equation |14[) 



rQ+ -— B* 
B-) 



(25) 



(26) 



(27) 



The adopted uniform distribution of B^ is consistent with 
transverse mechanical equilibrium at the inlet. We have also 
experimented with nonuniform distributions of the magnetic 
field, in particular with B^ decreasing with ^. The results 
were not significantly different as the field distribution down- 
stream of the inlet underwent a rapid rearrangement that 
restored the transverse force balance. 

To have control over the mass flux, the flow at the inlet 
boundary is set to be super-slow-magnetosonic. This means 
that both the density and the radial component of the ve- 
locity can be prescribed some fixed values: 

In the simulations we use Vpg — 0.5 c or 0.7 c, which is a 
choice of convenience. On the one hand, this value is suffi- 
ciently small to insure that the fiow at 77^ = 1 is sub-Alfvenic 
and hence that the Alfven and fast-magnetosonic critical 
surfaces are located downstream of the inlet boundary. On 
the other hand, it is large enough to promote the rapid es- 
tablishment of a steady state (in which the outfiow speed 
remains constant along the symmetry axis). Because of the 
sub-Alfvenic nature of the inlet flow, we cannot fix the other 
components of the magnetic field and the velocity — they 
are to be found as part of the global solution. Following the 
standard approach we extrapolate B"^ and B^ from the do- 
main into the inlet boundar y cel ls. We then compute v** and 



from equations ( 26 1 and ( 27 1 



The magnitude of the angular velocity is chosen in such 
a way that the Alfven surface is encountered close to the 
source. Specifically, in the case of solid-body rotation the 
light cylinder radius, ric, is — 50% larger than the initial jet 
radius. In the differential rotation case, the closest point of 
the Alfven surface is located at a distance of ~ 1 initial jet 
radius from the inlet surface. 

The inlet density varies from model to model in order 
to cover a wide range of initial magnetizations. Table I gives 
the key parameters of all the jet models constructed in this 
study. Most of the models, denoted by the letter B, corre- 



spond to the wall shape z oc 



p3/2 



and differ only by the 



value of the magnetization parameter: fj, varies from the rel- 
atively small value of 39, which is more suitable to AGN jets 
(see Paper I), all the way up to 620. Model B2H is included 
to study the effects of a high temperature at the source. 
The initial effective thermal Lorentz factor in this model 
is Fto = wo/poc? = 55. Models A and AW have a wall of 
conical shape. In model AW the half-opening angle of the 
cone is 90°, which allows us to model the case of an un- 
confined outflow (which could be relevant to pulsar winds). 
The remaining models help to explore the effects of differ- 
ential rotation (model D), of various other paraboloidal wall 
shapes (2 oc in model C, z ocr'^ in model F) and of a wall 
shape whose opening angle increases with distance (model 
E). 
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Figure 1. Computational errors for models A (top row) and Bl (bottom row). The plots show the flow parameters fe(*), f2(*) and 
/i(*) at the inlet (solid lines) and at ?? = 1 x 10^ for model A and 77 = 5 x lO'^ for model Bl (dashed lines). 



Table 1. Parameters of simulation models. 



Model 


a 


rotation 






or 6j 


^max 


A 


1 


uniform 


1.0 




= 0.2 


560 


AW 


1 


uniform 


1.0 




= 7r/2 


560 


Bl 


3/2 


uniform 


1.0 




= 2.0 


620 


B2 


3/2 


uniform 


1.0 




= 2.0 


310 


B2H 


3/2 


uniform 


55 




= 2.0 


370 


B3 


3/2 


uniform 


1.0 




= 2.0 


155 


B4 


3/2 


uniform 


1.0 




= 2.0 


78 


B5 


3/2 


uniform 


1.0 




= 2.0 


39 


C 


2 


uniform 


1.0 




= 2.0 


620 


D 


3/2 


(liHrrcntial 


1,1) 




= 2.0 


()()() 


E 


2/3 


uniform 


1.0 




= 0.1 


300 


F 


3.0 


uniform 


1.0 




= 2.0 


540 



3.1.2 Other boundaries 

The computational domain is always chosen to be long 
enough for the jet to be super-fast-magnetosonic when it 
approaches the outlet boundary rj = r)o. This justifies the 



use of radiative boundary conditions at this boundary (i.e. 
we determine the state variables of the boundary cells via 
extrapolation of the domain solution) . 

At the polar axis, ^ = 0, we impose symmetry boundary 
conditions for the dependent variables that are expected to 
pass through zero there, 

f{-i) = -fio ■ 

These variables include B^, B^, and u'^. For other vari- 
ables we impose a "zero second derivative" condition, 

d\f/de = o, 

which means that we use linear interpolation to calculate 
the values of these variables in the boundary cells. 

We do this in order to improve the numerical represen- 
tation of a narrow core that develops in all cases as a result 
of the magnetic hoop stress. Within this core the gradients 
in the ^ direction are very large and the usual zero-gradient 
condition, f{—^) = /(Oi results in increased numerical dif- 
fusion in this region. We have checked that this has a no- 
ticeable effect only on the axial region and that the global 
solution does not depend on which of these two conditions 
is used. 

At the wall boundary, ^ = ^j, we use a reflection condi- 
tion, 
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Figure 3. Same as in Fig.|2] but for model D. Tlie closest to the inlet point of the Alfven surface has the radius ri^ = 1.3. 



S. S. Komissarov et al. 




Figure 5. Same as in Fig. [2] but for model B2H. Tlie light cylinder radius is rjc = 1.6. 
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for and and a zero-gradient condition for all other 
variables. 



3.2 Initial setup 

The initial configuration corresponds to a non-rotating, 
purely poloidal magnetic field with approximately constant 
magnetic pressure across the funnel. The plasma density 
within the funnel is set to a small value so that the out- 
flow generated at the inlet boundary can easily sweep it 
away. In order to speed this process up the r) component of 
velocity inside the funnel is set equal to 0.7 c, whereas the ^ 
component is set equal to zero. 

3.3 Grid extensions 

The inner rings of the grid, where the grid cells are small 
and, therefore, so is also the time step, are the computation- 
ally most intensive regions of the simulation domain. If we 
kept computing these inner rings during the whole run then 
we would not be able to advance very far from the jet origin. 
Fortunately, the transonic nature of the jet flow allows us 
to cease computations in the inner region once the solution 
there settles to a steady state. To be more precise, we cut 
the funnel along the ^-coordinate surfaces into overlapping 
sectors with the intention of computing only within one sec- 
tor at any given time, starting with the sector closest to 
the inlet boundary. Once the solution in the "active" sector 
settles to a steady state we switch to the subsequent sec- 
tor, located further away from the inlet. During the switch 
the solution in the outermost cells of the active sector is 
copied into the corresponding inner boundary cells of the 
subsequent sector. During the computation within the lat- 
ter sector these inner boundary cells are not updated. This 
procedure is justified only when the flow in a given sector 
cannot communicate with the flow in the preceding sector 
through hyperbolic waves, and thus we ensure that the Mach 
cone of the fast-magnetosonic waves points outward at the 
sector interfaces (see Paper I). 

In these simulations we used up to 7 sectors, with each 
additional sector being ten times longer than the preceding 
one. This technique has enabled us to reduce the computa- 
tion time by more than three orders of magnitude. Although 
the grid extension can in principle be continued indefinitely, 
there are other factors that limit how far along the jet one 
can advance in practice. Firstly, once the paraboloidal jets 
become highly coUimated the required number of grid cells 
along the jet axis increases, and each successive sector be- 
comes more computationally expensive than the previous 
one. Secondly, errors due to numerical diffusion gradually 
accumulate in the downstream region of the flow and the 
solution becomes progressively less accurate (see Fig. [ij. 



4 RESULTS 

As is generally the case in numerical simulations, our com- 
putations are subject to numerical errors, mainly the trun- 
cation errors of our RMHD scheme. The field-line constants 
described in Section |2.1| can be used for a straightforward 
evaluation of the absolute error. Fig.[T]shows the ideal-MHD 
constants k, Q. and fj, as functions of magnetic flux at the 




0,5 



Figure 6. Distribution of the poloidal magnetic field across the 
jet of model Bl, showing the development of an axial core as the 
distance from the origin increases. From top to bottom, the curves 
correspond to ?? = 1, 50, 5 X 10^, 5 X 10^, 5 X lO", 5 X 10^, 5 x 10^ 
and 5 X 10^, respectively. 



inlet and near the outer boundary of the computational do- 
main for models A and Bl. If the curves do not exactly coin- 
cide, this is indicative of computational errors. Although the 
plots exhibit noticeable deviations, they remain relatively 
small, and we conclude that the results are trustworthy. 

Figs. [2}j5] show the general 2D structure of the derived 
jet solutions for models A, D, B2 and B2H midway from the 
inlet surface. We selected these particular cases since they 
represent the most significant variations in the model param- 
eters, namely the transition (i) from conical to paraboloidal 
shape of the confining wall (A and B2), (ii) from uniform to 
differential rotation at the base (A and D|^ and (iii) from 
cold to initially hot flows (B2 and B2II). In general, the 
structure of the simulated ultra-relativistic jets is very simi- 
lar to that of the moderately-relativistic conical jets studied 
in Paper I. All models show the development of a central 
core where the source-frame mass density p' — Tp peaks. 
The mass concentration is accompanied by a buriching-up 
of the poloidal magnetic field lines near the axis, as further 
illustrated in Fig. |6] The development of an axial core is a 
generic property of axisymmetric MHD outfiows from a ro- 
tating source ( Bogovalov|1995 \ and was also a feature of the 
jets simulated in Paper I. The distribution of the Lorentz fac- 
tor across the jet varies, however, from case to case. In model 
A r has its maximum value at the jet boundary (Fig. [2|. 
In model D the maximum is located approximately mid- 
way between the symmetry axis and the boundary (Fig. [3|. 
This reflects the fact that the angular velocity of magnetic 
field lines, and hence the electromagnetic energy flux (equa- 
tion 201, vanishes at the boundary in this model, resulting 



* Note that, when displaying results for model D, we define the 
fiducial light-cylinder radius in terms of the angular velocity Qq 
of the innermost field line. 
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Figure 7. Distribution of F and /irn, = MfeO" across the jet in models Bl (left panel) and D (right panel). Solid Unes show F at 
= 5 X 10*, 5 X 10^, 5 X 10®, 5 x 10^ (increasing upward), dashed Unes show fii^a at the same locations (increasing downward), and the 
dash-dotted line shows /i. 



in /X ~ /i/i 1 near the wall (see Fig. [7|. The Lorentz fac- 
tor of the initially cold jet in model B2 at first peaks near 
the axis, with its value decreasing slightly on the way to the 
jet boundary. However, further downstream the maximum 
shifts towards the boundary and eventually disappears. In 
the initially hot jet of model B2H the Lorentz factor at first 
peaks right on the symmetry axis, where the acceleration is 
due to the by the gas pressure. However, further downstream 
its evolution is similar to that of model B2. 

Fig. [8] shows the efficiency of plasma acceleration along 
the magnetic surface 5' — 0.8^'max (located near the jet 
boundary) for models B1-B4, which differ only by the 
strength of the initial magnetization. One can see that in 
all four cases the kinetic energy flux, ~ fihpUpC^ — TpUpC^, 
eventually exceeds the Poynting flux, fimpUpC^. This mag- 
netic surface is not exceptional and a similar behaviour is 
exhibited along other flux surfaces. This is illustrated by 
Figs. [T] and [9] These figures also show that soon after reach- 
ing equipartition the plasma acceleration slows down signif- 
icantly: this is consistent with the relation jj, ~ r(l -I- a) 
obtained from equations ( 16 1, ([20| and il9l, in which cross- 



ing the equipartition point corresponds to the magnetization 
parameter a dropping below 1. Fig. [S] further indicates that 
the efficiency of magnetic acceleration is higher the lower 
the initial magnetization. This is reflected in the behaviour 
of a, the ratio of the Poynting flux to the matter energy flux 



(see Section 2.11. The left panel of Fig. 10 shows that the 
fast initial decrease of a slows down at a higher value of a 
when the initial magnetization is larger. If this behaviour 
in fact extends to values of jj,mO ~ P that are low enough 
for the maximum attainable speed to remain nonrelativis- 
tic then the indicated inverse correlation is consistent with 
the very high acceleration efficiency exhibited by MHD out- 
flow solutions in the Newtonian regime (e.g. [Vlahakis et al.| 



The high efficiency of magnetic acceleration is not 
unique to models in which the magnetic field lines rotate 
uniformly. Fig. [7] in which the results for model Bl are com- 
pared with those for model D, shows that equally effective 
acceleration is achievable in the case of a differentially ro- 
tating source. 

The geometry of the bounding wall has a pronounced ef- 
fect on the acceleration efficiency, as demonstrated by Fig.|9] 
A larger value of the power-law index a in the shape function 
z oc corresponds to a more rapidly rising function r(r/ric) 
along a given magnetic flux surface "if =const. Whereas in 
the model Bl (a = 3/2) the acceleration slows down only 
after the equipartition point, in model A (a = 1) this occurs 
much earlier and, as a result, equipartition between mag- 
netic and kinetic energy is reached only near the jet axis. 
Equipartition is not reached in model C {a — 2) either (see 
Fig.[9|, but for a different reason. Due to the higher degree of 
external coUimation, this jet eventually becomes very thin. 
This makes our simulation increasingly expensive and we 
are forced to terminate it before reaching sufficiently large 
jet radii. (Moreover, the computational errors are accumu- 
lated over a longer path along the jet and would become 
rather high if we continued.) However, Fig. [9] shows that in 
this model the Lorentz factor is a faster growing function of 
cylindrical radius compared to model Bl. Finally, in model 
E {a — 2/3) we consider a jet propagating in a channel with 
a progressively diverging wall, which in practice may cor- 
respond to the polar funnel of a thick accretion disc (e.g. 
Paczyhsky fc Wiita|198"o| . In this case the jet eventually be- 



comes detached from the wall and then expands as a conical 



outflow (Fig. 111. The acceleration rate is similar to that of 



2000) 



model A (see Fig. 121. 

The initially hot jet, model B2H, is subject to both mag- 
netic and thermal acceleration, so, as expected, the Lorentz 
factor in this case grows faster compared to the correspond- 
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Figure 8. F (solid line), fim = Mh"" (dashed line) and fj, (dash-dotted line) along the magnetic field line with = O.S^max as a function 
of cylindrical radius for models Bl (top left panel), B2 (top right panel), B3 (bottom left panel) and B4 (bottom right panel). 



ing cold jet (see the right panel of Fig. 131. But a closer 



inspection reveals that the acceleration process exhibits a 
new mode of behaviour in this case (one that was, however, 
found before in semi-analytic self-similar solutions; see |Vla-| 
hakis fc Konigl||2003b| ). It is seen that a significant fraction 
of thermal energy is at first converted into Poynting flux. 
The middle panel of Fig. [13] shows that the Poynting-to- 
mass flux ratio /i,„c^ grows until r ~ lO^ric and only then 
starts to decline. However, this decrease is quite fast and the 
terminal value of fim for the chosen magnetic flux surface 
(^f = O.S^max) is, in fact, lower than in the corresponding 
cold jet (model B2) shown in the left panel of this figure, 
with a correspondingly higher asymptotic Lorentz factor. 

The distribution of the terminal bulk Lorentz factor 
across these two jet models is shown in right panel of Fig. |13| 
One can see that on the axis the Lorentz factor of the hot 



jet is higher than that of the cold jet by approximately the 
value of the initial thermal Lorentz factor, Fto = 55. This is 
as expected given that magnetic acceleration does not oper- 
ate along the axis. However, at the wall the difference is only 
half as large and in the middle of the jet it is higher than 
40. These traits are evidently a consequence of the thermal- 
to-Poynting energy conversion and its effect on the poloidal 
magnetic field distribution, as discussed in Section [5. 5| 



Although the case of an unconfined wind may not be di- 
rectly relevant to GRB fiows, which are inferred to undergo 
a fairly efficient collimation (see Section [T]), it is certainly of 
interest to the pulsar community. Furthermore, it is worth 
investigating from a purely theoretical point of view. The ac- 
celeration details for this case (model AW) are presented in 
Fig. |14[ The lower efficiency of magnetic acceleration noted 
in the conical- wall case (model A), particularly near the jet 




—I I ; \ ; \ I O I 1 ' 1 ' 1 ' 1 
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Figure 10. Left panel: Evolution of a along the magnetic field line VP = 0.8<I/niax in models Bl (solid line), B2 (dashed line), B3 
(dash-dotted line), B4 (dotted line) and B5 (dash-triple-dotted line). Right panel; Evolution of the bunching function S = nBpr^ /'i for 
the same models along the same magnetic field line. 



boundary, is even more pronounced in this instance. As can 
be seen in the right panel of Fig. |14[ only ~ 5% of the Poynt- 
ing flux injected at ~ 12° to the equatorial direction has 
been converted into kinetic energy by the time the cylindri- 
cal radius grew to r = lO^ric. Although, as shown in the left 
panel of Fig. |14[ the efficiency is higher near the symmetry 
axis, the terminal Lorentz factor there remains compara- 
tively low because of the reduced effectiveness of magnetic 
acceleration as the polar angle approaches zero. 



5 ANALYSIS OF THE RESULTS 

5.1 Efficiency of magnetic acceleration 

The steady-state structure of a magnetized relativistic out- 
flow can be understood by analysing the momentum equa- 
tion. After the partial integration described in Section |2.1[ 
two more equations remain to be considered, corresponding 
to the two components of the momentum equation in the 
poloidal plane. Since the main part of the acceleration oc- 
curs in the super-Alfvenic region of the flow, it is sufficient 
to examine only this regime. We further simplify the dis- 
cussion by taking the flow to be cold. Thermal effects, when 
present, in any case only affect the initial acceleration region 
of the flow; we consider them in Section fS.SI We now proceed 
to extend the discussion in Paper I by taking the F ^ 1 of 
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the constituent equations, appropriate for the ultrarelativis- 
tic flows simulated in the present work, which enables us to 
derive analytic scalings. 

For cold flows /i^ ~ F (equation 19 1, and from equa- 
tion ( 16 1 one finds that F ~ /i — /im- Substituting the electric 



current from equation (241 into equation (20 1, we get 



where 

/ Bp-dS 



Tvr^Bp _ r|V*| 



2* 



Thus, the flow Lorentz factor can be written as 



o 5. 



(28) 



(29) 



(30) 



All the quantities except for S on the right-hand side of this 
equation are field-line constants, so an increase in F along 
a field line necessarily requires 5 to decrease. The function 
5 is a measure of how bunched the poloidal field lines are 
— indeed, it is equal to the ratio of Bp at some cylindrical 
radius r along the field line to the mean magnetic field within 
that radius, ^'/vrr^. For example, for a fiow confined within 
a sufficiently small angle that satisfies Bp cx r"^, *!/ cx r"^^^ 
and 



5 = 



A + 2 



For a uniform distribution of Bp this yields 5 = 1, whereas 
one has 5 > 1 if -Bp increases with r and 5 < 1 if it de- 
creases. This shows that magnetic acceleration requires a 
gradual concentration of magnetic flux in the central part 
of the flow. In the case of a collimating flow this can be 
achieved through a faster coUimation of the inner magnetic 
flux surfaces than of the outer ones, and in the case of a 
decollimating flow a faster decollimation of the outer flux 
surfaces is required. Fig. |6] illustrates the concentration of 
magnetic flux toward the axis in one of our simulations. In 
this case, at large distances the poloidal magnetic fleld scales 
roughly as Bp oc r~^'^, corresponding to 5oo ~ 0.4. This is 
indeed the asymptotic value of 5, as shown in Fig. |15| 

Equation ( |30[ ) is a consequence of the momentum equa- 
tion along the flow. It shows how F increases by the action 
of the (l/c)Jp X force when the function 5 decreases 
along the flow, thereby demonstrating the intimate connec- 
tion between the acceleration efficiency and the evolution of 
the poloidal shape of the flow. In evaluating this efficiency 
we can use Si, the value of S at the fast surface, as a con- 
venient proxy for the initial value of 5. This is because, for 
/i ^ 1, F remains ^ /x on this surface (e.g. Komissarov 



|2004[ ). In this case the two terms on the right-hand side of 
equation (30 1 are comparable, and we obtain 



St = 



We can legitimately use equation ( 30 1 since the fast surface 



lies well outside the light cylinder and hence is in the super- 
Alfvenic domain for most of the simulated field lines. We 
now utilize this equation to write the asymptotic Lorentz 
factor in the form 




-15 



-10 



Figure 11. Colour image shows logjQ ptot (with the total pressure 
given by ptot = P + ^co/^''"' where Bco is the eomoving magnetic 
field) and the contours show the magnetic field lines for model E. 
In this model the light cylinder radius is ri^ = 0.29. 




Figure 12. Lorentz factor along the magnetic field line with 
^ = O.S^niax for model A (solid line) and model E (dashed line). 



' m(1 ~ Soc/Si) 



(31) 
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Figure 13. Effects of thermal acceleration. Left panel: cold jet of model B2. Middle panel: hot jet of model B2H (with uiq/pqc^ = 55). 
The lines show F (solid line), fi (dash-dotted line), = fif^a (dashed line) and {wj p(? — l)r (dotted line) along the magnetic field line 
with 'I' = 0.5^max as a function of cylindrical radius. Right panel: T r^r.or^t., fc^t.^^ ar-moo +Vio lot c+ »i — /i v in6^ 
model B2 (solid line) and the hot jet of model B2H (dashed line). 



Lorentz factor across the jet at = 4 X 10 for the cold jet of 




Figure 14. Unconfincd wind solution (model AW). Left panel: Lorentz factor (increasing function) and /Jfecr (decreasing function) along 
five different magnetic field fines: * = O.S'I'max (sofid line), >!' = 0.5'I'max (dashed fine), * = 0.2VI'niax (dash-dotted fine), vl/ = 0.1'I'max 
(dotted line), f = 0.027'!' max (dash-triple-dotted line), the last line originating from the same point at the inlet as the ^ — O.S'I^niax 
line of model A. Right panel: F (solid line), ^;,o" (dashed line) and (dash-dotted line) along the magnetic field line with 'I' = 0.8'I'max 
as a function of cylindrical radius. 



In our simulations Si ~ 0.9 (see Figs. 10 and 15 1. This value 
reflects the adopted uniform distribution of _B'' at the in- 
letj^Beyond the Alfven surface the azimuthal magnetic field 
component becomes dominant, and its hoop stress causes 



^ As we already in Section [3.1.l| we have experimented with other 
distributions that put more flux near the axis and observed a 
quick "uniformization" of magnetic fiux in the immediate vicinity 
of the inlet under the action of magnetic pressure. 



the inner flux surfaces to collimate faster than the outer 
ones. As a result S decreases, attaining asymptotic values 



0.25 — 0.4 for paraboloidal jets (see Figs. 10 and 15 1 



The implied asymptotic Lorentz factors thus satisfy 
r^//i^ 0.55 - 0.72, 

which are indeed the values reached by our simulated flows 
(see Figs.[7]-[9|. This result indicates that > 50% of the initial 
Poynting mix is converted into kinetic energy of bulk motion 
(see also |Vlahakis|2004b I. The signiflcantly lower efficiency 
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Figure 15. Evolution of the function S = irBpr^ along the 
magnetic field line with = 0.5'!' max in models A (solid line), 
Bl (dashed line), C (dash-dotted line) and D(dotted line). 



found in our simulations of flows inside conical and diverging 
funnels, down to 25% near the boundary (models A and E), 
is most likely due to the loss of causal connection across the 
flow (see Section [5. 4[). 



5.2 Power-law acceleration phase 

Next we analyse the trans-fleld component of the momentum 
equation. The asymptotic form of the trans-field equation in 
the highly relativistic limit is 



7^ 





^ rVln 
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V* 




f 






1 + 


W 4:TVpU 


p T 


2 

Ic 






T 


2 





-.2n< 



Vr-V* 



|V*| 



(32) 



where TZ is the curvature radius of poloidal field lines (see 
equation 16 and related discussion in Vlahakis ''2004a ) . The 
three terms of this equation are the poloidal curvature term 
(left-hand side), the electromagnetic term (first on the right- 
hand side), which is of order 1, and the centrifugal term (sec- 
ond on the right-hand side). This important equation, with 
the centrifugal term omitted, was derived by|Chiueh et al.] 
|1991l),|Lyubarsky & Eichlerl (120011), andlOkamotol (120021) , 



while |Bogovalov| (|1995[), [Beskin fc Malyshkin| (|2000|, and 
[Tomimatsu fc Takahashi| ( |2003[ ) derived the same equation 
with the centrifugal term included but the poloidal curva- 
ture term omitted. 

Well outside the light cylinder, where rQ, ^ u*^ and 
V ~ c, equations ( 26 1 and ( 27 1 imply 



rS'' = --nBpr^ . 
c 



From this equation and equation ( 18 1 one finds that 
1. 



-QB^r 



(33) 



(34) 



where Bp is the magnitude of the poloidal magnetic field. 



Substituting this result into equation ( 20 1 one also finds that 
1 



b;v 



47r rf^ pUp 



Thus, in this regime one can rewrite equation (32 1 as 
f 



rVln 
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1 + 



27^ Vr-V* 
r2 |V*1 



(35) 



(36) 



In the magnetically dominated case, where p.m 2> p-h, order- 
of-magnitude evaluation of the last two terms in this equa- 
tion gives the useful result 



71 



1 



p2'i 



(37) 



Depending on which term in equation ( |32[ ) can be ne- 
glected, we can isolate the following three cases (ordered by 
increasing importance): 

(i) If the electromagnetic part is negligible then the 
shape of the flow is determined by the centrifugal term, re- 
sulting in a hyperbolic line shape, a characteristic of ballistic 
motion (see equation 20 and related discussion in |Vlahakis| 
2004a see also Sections 5.3.1 and A1.3I. None of the end- 



states of our simulations has this property. 

(ii) If the poloidal curvature term is negligible, the elec- 
tromagnetic and centrifugal terms balance each other. This 
is the case very close to the rotation axis (inside the cylin- 
drical core) as well as for a quasi-conical flow like our model 
A and for paraboloidal fiows with a > 2 as in our model F 



(see Section 5.3). In this case equation (37 1 gives 



(38) 



Following different methods, this "linear acceleration case" 
was found by [Contopoulos fc Kazanas] ( |2002[ ), who anal- 
ysed radial force-free flows beyond the light cylinder (and 
hence their analysis holds in the regime between the Alfven 
and the fast-magnetosonic surfaces), and by [Beskin et al.] 



( 1998 1, who perturbed a quasi-conical flow (and found that 
F « r/ric applies in the sub-fast-magnetosonic regime). Our 
results for models A and F agree with the scaling F ~ r/ri^; 
see the top left panel of Fig. |17| 

(iii) If the centrifugal term is negligible then the shape 
of the flow is determined by the electromagnetic force. This 
regime applies to the case of paraboloidal wall with a < 2 



(see Section 5.3 1. Equation (37 1 implies that in this case the 



(39) 
(In what 



radius of curvature of poloidal field lines is 
7^^FV. 

Now, consider a field line of the shape, z a 
follows we use the superscript b to indicate the power-law 
index that describes the shape of given magnetic field lines, 
whereas the superscript a is reserved for the power-law index 
that gives the shape of the funnel wall in our numerical 
models. Note that the interior field lines in these models 
have b that is slightly larger than a, although & — + a as the 
wall is approached; see Fig. 16 ) The curvature radius of 
such a line satisfies 



r 
7^ 



a r 



&2 



(40) 
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where the final form is valid when Bp ^ Bz. Combining this 



with equation ( 39 1 we get 



Z 6-1 (6-l)/6 



(41) 



(see also Vlahakis fc Konigl|2003b I, which applies when the 
power-law index lies in the range 1 < 6 < 2 and shows that 
the spatial growth of the Lorentz factor is also a power law 
in this case (in either r or z). Assuming that the flow is not 
too collimated within the light cylinder, so that zic — ric for 
most of the fleld lines (an assumption that is well satisfied 
in our numerical models), we can write the above result in 
the following useful forms: 



(r/ric 



x(6-l)/6 



(42) 



This acceleration regime operates in our 1 < a < 2 numerical 
models before the fiow reaches approximate equipartition, as 
can be verified by inspecting Figs. |17|and |18| 

The direct dependence of the flow acceleration on the 
poloidal curvature of the magnetic field lines in the regime 
(iii) leads to an anti-correlation between the jet Lorentz fac- 
tor and its opening angle. For a line shape 2:ocr''(l<6<2) 
we find 



rtan^v = 1/Vb- 1 , 



(43) 



where 6^ = arctan(dr/ci2) is the local half-opening angle 
of the magnetic flux surface. Fig. [l9] shows the variation of 
rtan^v along the flux surface — CS^tmax of model Bl. 
One can see that this product is indeed close to 1/ \/& — 1- It 
is, however, not exactly a constant, for the following reasons: 
the curvature acceleration regime is not really applicable 
at small and large spherical radii, the electromagnetic term 



in equation (321 is not exactly equal to 1, and the power- 
law index b varies along the fiow. The figure nevertheless 



indicates that equation ( 43 1 provides a useful estimate of 



the relationship between F and 9^. 

As expected from our discussion in Section |5.1| of the 
close connection between the acceleration efficiency and the 
evolution of the poloidal field-line shape, the trans-field force 
balance equation, which determines the variation of the 
flux-surface shape along the flow, is seen to provide infor- 
mation on how fast the Lorentz factor increases with dis- 
tance from the source. For all shape functions z x r'' with 
1 < 6 < 2, the corresponding power-law dependence of F 
leads to a high ( > 50%) magnetic-to-kinetic energy conver- 
sion efficiency over astrophysically relevant distances. Using 
equation ( 42 1 we find that equipartition between the Poynt- 



ing and kinetic energy fluxes is attained at a cylindrical ra- 
dius 



ro 



(2F0) 



1/(6-1) 



(44) 



After substitution ro = ric and Fq = 1, this equation reads 

/U\ 1/(6-1) 

^c, = nc(|j , (45) 

which, in fact, agrees very well with our results for models 



B (see Fig. 191. In terms of the spherical radius, assuming 



again that Ric — ric, we can write this expression as 

^e, = nc (f ) . (46) 



For b > 2 the corresponding relations are (using equation 38 1 
rcq = (p/2)ric and 7?cq = {fi/2fric. 

The derived scaling for the Lorentz factor can be used to 
flnd the behaviour of other quantities. For example, for the 
main part of the flow in which the Foynting flux dominates 
the energy flux, one has crF ~ ^ and hence, for 1 < 6 < 2, 



n/T oc r/ 



oc r z oc r 



-(6-1) 



(47) 



The predicted behaviour is indeed seen in the left panel of 
Fig. |10| This figure further shows that the "self similar" 
structure of the magnetization curves extends also beyond 
the equipartition radius, where they flatten out; in particu- 
lar, they do not cross each other even in that regime. Con- 
sequently, the magnetization beyond the turning point of 
the curve is lower the smaller the inlet value, which goes 
along with our flnding that the efficiency ~ 1/(1 + a 00) of 
magnetic-to-kinetic energy conversion in cold flows decreases 
with increasing initial magnetization. 

The high acceleration efficiencies attained by our sim- 
ulated flows appear to be inconsistent with the conclusion 
of |Chiueh et al.| ( |1998| ) that a transition to a low-cr conflg- 
uration cannot occur gradually in regions well beyond the 
light cylinder, where the fiow has become ultra-relativistic. 
Their analysis was, however, based in part on an estimate 
of the change in the angle 9^ between the poloidal flow and 
the rotation axis as one moves through a length along 
the fiow (see text after equation 14 in their paper): this es- 
timate is not generally valid since it assumes that ~ Ar, 
which only applies to quasi-radial flows. If instead we use 
Al ~ Az in equation (14) of Chiueh et al. (19981 and con- 
centrate on paraboloidal flows (z oc r ) with & < 2, we get 
ASv ~ A(z/r)6/F^(fo - 1), which yields the scaling F oc z/r 
found above. On the other hand, the lower acceleration ef- 
flciency exhibited by our model A, in which the flow mor- 
phology is quasi-rad ial (see Figs. [2} [9| and 15 1, appears to be 
" " " " |l998l " 



consistent with the Chiueh et al. (11998 1 inference of loga- 



rithmic coUimation and slower acceleration. We note in this 
connection that, beyond the end of the power-law accelera- 
tion phase analysed in this subsection, it is possible to have 
an additional, logarithmic acceleration regime in which po- 
tentially up to 100% of the Foynting flux could be converted 
into matter kinetic energy flux (see Vlahakis'2004a' and ref- 
erences therein). However, this acceleration is too slow to be 
of astrophysical interest since it requires exponentially large 
distances for completion. 



5.3 Dependence on the external pressure 
distribution 

Although we have chosen, for numerical convenience, to pre- 
scribe the shape of the funnels that guide our simulated 
flows, in reality the boundary shape of pressure-conflned 
flows will be determined by the ambient pressure distribu- 
tion, pcxt, and we expect a one-to-one correspondence be- 
tween the shape of the boundary and the parameters of the 
confining medium, enforced through the pressure-balance 
condition at the boundary, pint = Poxt ■ Here we analyse this 
issue for the asymptotic region of a magnetically accelerated 
flow, where the internal jet pressure, pint, is dominated by 
the contribution due to the azimuthal component of mag- 



Magnetic acce. 



leration of ultra- relativistic jets 

^ ro I — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — , — ,-, 



in gamma-ray burst sources 

kD I n 1 1- 1 1 1 1 1 1 1 1 1 1 1 1 



17 



Figure 16. The exponent b of the poloidal shape function z oc for models A (left panel), B2 (middle panel) and B2H (right panel) 
across the jet. For model A the depicted cross sections are at R = 10 (solid line), _R = 10^ (dashed line), R = 10'^ (dash-dotted line), 
R = 10* (dotted line) and R = 10^ (dash-triple-dotted line). For models B2 and B2H the plotted cross sections are at r; = 5 X lO'^ (thin 
solid line), t; = 5 X lO'^ (dashed line), r/ = 5 X 10** (dash-dotted line), rj = 5 X 10^ (dotted line), r] = 5 X 10® (dash-triple-dotted line) and 
= 5 X lO'^ (thick solid line). 



netic field, pint =p + B^^/8n ~ (B'*)^/87rr^ Thus, 

^ (B<^)2 ' 

In the following we assume that the external pressure dis- 
tribution is a power-law 

Poxt = Pc,it,lc{z/z\cy" , 

which is consistent with the funnel shape 2; oc adopted in 
our numerical simulations. Moreover, since fj,m oc / oc rB'^ 
(see equations |18[[20| ) is a weak function of distance we may 

assume that at the jet boundary B''' = Bf'^{r /ric)"^ , Then 
we have 



(48) 



where x = r/ric and Z = z/ric are the dimensionless coor- 
dinates of the jet boundary and 



C = 



(zic/ricY 

*- Ic 



(49) 



It is easy to see that C is a positive dimensionless constant of 
the order of 1. Provided that dr/dz ^ 1 we can approximate 
the curvature radius of the jet boundary via 



d^r 

'd^ 



1 d^x 
ncdZ^ 



and rewrite equation (37 1 as 
d'^x 



' dZ^ 



1 

f2 



1 



0. 



(50) 



(51) 



After the substitution of F from equation ( 48 1 this yields an 



ordinary differential equation for the jet boundary 



d X X 
dZ^ Z'- 



= 0. 



(52) 



The first term on the left-hand side of equation ( 52 \ rep- 



resents the effect of poloidal curvature, the second is the 
electromagnetic term and the third is the centrifugal term. 



Equation ( 52 \ can be solved in closed form in various 



limits, as discribed in Appendix |X] Here we simplify the 
discussion by looking for almost power-law solutions 

x^K-^Z^^", (53) 

with K being positive constants and a varying very slowly. 



Substituting this ansatz into equation ( 52 1 and ignoring all 



terms including derivatives of a, we obtain 



-(--0 
a \a J 



CZ^ 



0. 



(54) 



We now proceed to analyse this equation for different values 
of the exponent a. 

5.3.1 a>2 

In this case the second term on the left-hand side of equa- 
tion ( 54 1 vanishes aa Z ^ 00 and the only acceptable asymp- 



totic value of a is unity. Indeed, for a > 2 the third term 
diverges, for a = 2 it is constant but negative and so is the 
first term, for o. < 2 it vanishes and so must the first one, 
implying a — *■ 1. Thus, asymptotically the boundary adopts 
conical shape. 



• When a < 4 the electromagnetic term of equation ( |54[ ) 
dominates over the centrifugal term, and thus a ^ T^ (since 
the first term must be negative in order to cancel the sec- 
ond). The boundary shape is therefore paraboloidal (with 



conical asymptotes). An explicit solution of equation ( 52 1 in 
this limit is given in Appendix [A] 

• When Q > 4 the centrifugal term dominates over the 
electromagnetic term in equation (54 1 and thus a ^ 1~ 



(since the first term must be positive in order to cancel the 
third). This is case (i) of our analysis of equation | |36[ ), which 
corresponds to a hyperboloidal shape (with conical asymp- 
totes), as demonstrated in Appendix [A| through an explicit 
solution of equation ( |52[ ) in this limit. 

• When a — 4 one can obtain a solution that is conical 
(a — 1) from the start, with Ji"' = C. This solution corre- 
sponds to our conical model A during the acceleration phase, 



when r oc r (see equation 381. Fig. 20 verifies the predicted 
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Figure 18. Lorentz factor along three different magnetic field lines of models A (top left panel), Bl (top right panel), C (middle left 
panel), D (middle right panel), F (bottom left panel), and B2H (bottom right panel) as a function of the spherical radius R. Solid line: 
* = 0.8*max; dashed line; * = 0.5*max; dash-dotted hne: * = 0.2*niax. 
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Figure 19. Left panel: variation of FtanSv along the flux surface 'I' = O.S'I'niax of model Bl. Right panel: the diamonds show the 
equipartition radius (where the Poynting and kinetic energy fluxes are equal) along = 0.8'I'max as a function of the magnetization 
parameter fi for models Bl— B5. The solid line shows the function logjg(r/ri(.) = 2 logjg(/^/2). 



scaling (pcxt oc Z"**) and also shows that, after the growth 
of r saturates, a conical shape can be maintained only if the 
ambient pressure scales as z~^, which follows directly from 



the scaling pin 
subsection. 



oc r r discussed at the beginning of this 



In summary, for a > 2 the boundary does not simply ad- 
just to the ambient pressure profile but instead asymptotes 
to a conical shape. This result is consistent with the expecta- 
tion that in this case the transverse expansion time of the jet 
becomes shorter than the propagation time of magnetosonic 
waves across the flow, leading to a loss of causal connectivity 



and hence to a "free" ballistic expansion in a cone (Begel- 
man et al.|1984 see also Section 5.41. This is essentially the 
behaviour exhibited by our Model E (see Fig. 



5.3.2 a = 2 

In this case the second term on the left-hand side of equa- 



tion (54 1 is a positive constant. This implies that 1 < a < 2. 
(Indeed, for a > 2, the third term diverges and hence un- 
balanced. For a < 1 it vanishes but the first term is non- 
negative and hence cannot balance the second one.) We can 
distinguish between the following two cases: 

• a = 2 — the power law solution with A'* = C — 1/4 is 
exact. This implies C > 1/4. 

• 1 < a < 2 — the third term becomes negligible at 
large Z and balancing of the first two terms requires a 
2/(1 + VI - 4C). This implies C < 1/4. 

In other words, for C < 1/4 the centrifugal term is neg- 



ligible and the resulting shape is Z — {zic/ric)x' 



2/{l + v/l-4C) 

whereas for C > 1/4 the centrifugal term is comparable to 
the other two terms and the solution is Z = \/ C ~ 1/4 . 
Fig. |20| verifies that the confining pressure in our simulated 
flows scales as Z^^ irrespective of the precise value of a so 
long as the shape exponent lies in the range 1 < a < 2. The 
flgure also corroborates the prediction that the Z~^ scaling 
is attained only gradually when a < 2 (models B and D, cor- 



responding to a — 3/2) but that it is present almost from 
the start when a — 2 (model C). 

As shown in Appendix |Aj the asymptotic solution for 
C = 1/4 is a; = Z^^^iCi + C2 InZ), where CI and C2 / 
are constants. (We kept the constant Ci to accommodate 
the possibility that the solution extends all the way down 
to the light-cylinder radius, where Z ^ 1.) This solution is 
similar to the C < 1/4 solutions of equation ( |52[ ) in having 
a negligible centrifugal contribution. 

Although all the funnel shapes whose power-law indices 
lie in the range 1 < a < 2 correspond to a single expo- 
nent (a — 2) of the conflning pressure distribution, there 
is nevertheless a one-to-one match between a given pres- 
sure distribution and the resultant funnel shape. This is be- 
cause both the power-law index a and the magnitude of the 
confining pressure (as expressed in relation to the internal 
magnetic pressure at the light-cylinder radius by the pa- 
rameter C; see equation 49 1 play a role in determining the 



functional form of the boundary: when C < 1/4 the mag- 
nitude of C flxes the exponent of the boundary paraboloid, 
whereas when C > 1/4 it fixes the normalization constant 
K. The parameter C is evaluated at the effective base of the 
asymptotic region of the flow and it conveys physical prop- 
erties (e.g. zic and Tie; see equation 49 1 imprinted on the 
outflow before it reaches the asymptotic regime. Thus, the 
asymptotic shape of a jet propagating through a power-law 
pressure distribution is determined both by the exponent of 
that distribution and by the evolution of the outflow before 
entering the asymptotic region. 



5.3.3 a<2 

In this case the second term on the left-hand side of equa- 



tion (541 diverges as Z — > 00. To balance this term, the 



third term must also diverge in this limit, which implies 
that a = 4/q > 2 and C = K'^. Thus, the jet shape is 
paraboloidal, Z = C^'^^a;*''". Like in the a = 2 case, both 
the parameters a and C are needed to uniquely flx the func- 
tional form of the jet shape. For q = 4/3 we have a — 3, the 
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funnel shape index of our numerical model F. Fig. |20| verifies 
that the boundary pressure for this model indeed scales as 

^-4/3^ 

We can collect the results derived in this subsection 
into a concise description of the correspondence between the 
exponent a of the ambient pressure distribution and the 
exponent a of the asymptotic jet shape: 

• a < 2 <4> a = 4/q > 2 , 

• a = 2 <4> l<a<2, 

• a>2 <4> a=i. 

Similar results for the behaviour of the ambient pressure in 



a confined jet {a < 2) were found by Tchekhovskoy et al. 
(20081 in the force-free limit, which is consistent with the 



fact that our expressions for the spatial profile of F were 
obtained in effectively the same approximation. 

As we have seen, a — 2 leads to the asymptotic balance 
between the electromagnetic and poloidal curvature forces 
(regime iii) whereas a < 2 leads to the balance between the 
electromagnetic and centrifugal forces (regime ii; see Sec- 
These regimes are characterized by different evo- 



5.21 



tion 

lution of many flow parameters, which may have observable 
consequences (see also Section [gJ. For example, in regime 
(ii) the product Ftan^v is predicted to be a constant 0(1) 
in the acceleration region, whereas in regime (iii) it is ex- 
pected to decrease with distance as Z~^^~'^^''\ with b being 
slightly larger than a due to the stronger collimation of the 
flow inside the jet. The evolution of the Lorentz factor in 
regime (ii) is given by F oc r (equation 38 1 rather than by 
the F oc r''~^ scaling of regime (iii). However, in practice 
this may not translate into a significant difference in how 
fast the jet accelerates (for example, F « z^^^ for both the 
a — 4/3 and a — 2, b — 3/2 cases). 

After the end of the acceleration the internal pressure 
scales as (since F = Foo = const). If the external pres- 
sure continues to decline as z~" , the pressure balance implies 
that the radial coordinate r increases faster compared to 
its variation during the acceleration. The new flow shape is 
Z — C^^"F^°a;^^" as a result of equation ( 48 1. For example, 
in the cases a = 2, 1 < a < 2, the flow becomes radial and 
the opening angle of the jet remains constant. The quantity 
F tan 9v is also constant and equal to C~^^^ = a/\/a — 1 (us- 
ing the relation between C and a, see Section [53^ . Thus, 
Ftan^v is a times larger compa red to its value during the 
acceleration phase (see equation 43 1 ^ 



5.4 Magnetic acceleration and causality 

We have found that the acceleration efficiency is smaller 
when the wall has a conical shape (model A) than in the 
cases when its shape is paraboloidal (see Fig. |9|. In the 
conical-wall case the flow attains equipartition only along 
field lines that are close to the rotation axis < 0.2^^ max). 



^ The change of this quantity is smooth and happens as the 
funct ion T{Z) changes from a power law to a constant. Equa- 
tion ijisjl, written as x = C-l/^Z"/^ [r(Z)]"\ gives Ftanfiv = 
rdx/dZ = C-i/2(q,/2 - dlnr/dInZ)Z°/2-i. In the cases with 
o = 2, l<a<2 the slope dlnV /din Z changes from 1 — 1/a 
during the main part of the acceleration phase (see equation |41| 
to zero after it ends. As a result, Ftanfiv changes from l/y/a — 1 
to a/\/a — 1. 
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Figure 20. Evolution of total pressure along the jet boundary in 
models A (solid line), Bl (dashed line), C (dash-dotted line), D 
(dotted line) and F (dash-double-dotted line). 



In accordance with our discussion in Section |5.1[ the vari- 
ation in the acceleration efficiency is tied to the difference 
in the degree of the collimation across the outfiow, as seen 
in Fig. |16| Only for small values of 5' does the exponent b 
become significantly larger than 1, corresponding to the in- 
nermost field lines bending toward the rotation axis, which 
implies that the bunching function 5 decreases along this 
portion of the outflow. In order for collimation to occur, 
there must, however, exist causal connectivity across the 
outflow. A related discussion of this issue can be found in 
Zakamska et al. (20081. However, the simpler flow structure 



assumed in that paper excludes the possibility of magnetic 
acceleration. In particular, the assumption of zero azimuthal 
speed implies that the current / is a constant of motion (see 
equation 15 I, which in turn means that firn remains constant 



(see equation 20 1. 

One can check whether the condition of causal connec- 
tivity is satisfied by comparing the field-line opening angle 



Bv (defined in Section 5.21 with the half-angle of the Mach 
cone of fast waves, 6n^. The latter can be found from the 
relation 



sin 9n 



FfCf 
Fw„ 



(55) 



where Cf and Ff are the fast speed and the corresponding 
Lorentz factor, respectively. Since FfCf — i3co/\/47rp, where 
Bco is the magnetic field as measured in the fiuid frame, and 
Vn ~ c, we have 



/ ^co 

I 47r pc^ 



1/2 



(56) 



In the magnetically dominated regime a ~ m/F. For highly 
super-magnetosonic fiows ^ 1. Thus, we may write 



^m-T^. (57) 

In the hydrodynamic limit the fast magnetsonic speed re- 
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duces to the sound speed and Ff , Cf in equation ( 55 1 should 
be replaced by Fs, Cg. For the ultra-relativistic equation of 
state and F ^ 1 this gives 8m — 1/F, the value used for 
causality analysis in |Zakamska et al. (20081. However, in 
the magnetic case 9m can be much higher because the mag- 
netosonic speed can be much closer to the speed of light. 
In the conical case we have F « R/ric, and 



(58) 



grows rapidly to a value > 1 (where it is a good approxi- 
mation to replace 9^ by 9). The left panel of Fig. 21 shows 
that only the inner part of the jet has 0v/^m < 1 and thus 
in causal connection. CoUimation (and thus efficient accel- 
eration) is possible only in this inner region. In contrast, the 
outer parts of the conical jet lack causal connection with the 
axial region and the flow there is essentially ballistic. 

In the paraboloidal case with b < 2 (for which ~ 1/F 
and F ^ {R/ri^Y''-^^^'') 



%/9m ~ (F/m)'/' ~ (l/M'/')(i?/nc)<''-'^/^"'' , 



(59) 



so this ratio grows much slower compared to the conical 
case. Moreover, the loss of causal contact formally occurs 
when F ~ /i, i.e. at the end of the acceleration phase. This 
is confirmed by our simulations. As one can see in the mid- 
dle and right panels of Fig. |21| during the power-law accel- 
eration phase 9v/9m grows slowly but remains less than 1 
almost everywhere in our numerical models. It subsequently 
decreases again when the growth rate of F goes down. 

In contrast, in the paraboloidal case with b > 2 (for 
which 9v « r/bz and F ~ r/r-ic), 



%/9n 



(l/6/''C"''*)(r/nc) 

(l/6M'^'c"/')(i?/nc) 



{5/2)-6 
(5/2b)-l 



(60) 



(see Section 5.3.31, and this ratio actually decreases with 
distance for b > 5/2! One can also argue quite generally 
that, even if F were to increase all the way up to /i, the value 
of the above ratio in that region, which can be estimated to 
be ~ l/bfi''~'^C''^'^, would likely remain < 1 (since 6 > 2, 
^ > 1 and C is of the order of 1; see equation 49 1. Thus, 



the necessary (but not sufficient) conditon for acceleration 
is satisfied in this case. This suggests that the acceleration 
efficiency may be comparable to the 1 < 6 < 2 cases. 

The behaviour of an unconfined wind is similar to that 
of an outflow in a conical funnel, which is not surprising 
given the fact that the former is a limiting case of the latter. 

the acceleration in model AW is > 50% 
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As seen in Fig. 
efficient only along field lines that are close to the rotation 
axis (vl> < O.l^&max), similarly to the situation in model A. 



5.5 Hot flovirs 

When w/pc? is significantly larger than 1 at the inlet there 
is an additional reservoir of energy for the flow acceleration 
— the thermal energy of particles. As the flow expands the 
enthalpy per unit rest mass w/p — c? + [■^/(•s — l)](p/p) 
(equation 111 decreases until it reaches its minimum value 



(= c ), and beyond that point the flow can be regarded as 
cold. In the pure hydrodynamic case the thermal energy is 
directly transferred to the bulk kinetic energy of the fluid. 
In the magnetic case there is an additional possibility — the 
thermal energy can also be transferred to the Poynting flux. 



Indeed, since p,(? = {w/p)V + p.m(? , it is possible to have 
both F and ^m(? increasing when w/p decreases, and this 



in fact is what we observe in model B2H (Fig. 13 1. 

We have already noted in Section |5.1| that in Poynt- 
ing flux-dominated flows is proportional to the bunching 



function S (see equation 28 1. In agreement with this result, 
the left panel of Fig. |22| shows that in model B2H S exhibits 
the same evolution as jim (which is shown in the middle 



panel of Fig. 13 1 



In the super-Alfvenic regime the trans-fleld force bal- 
ance for hot flows is described by equation ( 32 1 even for 
hot outflows provided that p remains ^ B^o/^tt. Therefore 
we still have TZ ~ F^r and hence F oc r^~^ along magnetic 
fleld lines of paraboloidal jets with exponents in the range 
1 < fe < 2. Combining the mass conservation relation ( 13 1 
and equation p9| we obtain 



Fp = oc r , 

where we took account of the fact that 5 is a weak function 
of r. This enables us to write the variation of the thermo- 
dynamic parameters as 



p oc r 



p (X r 



-s(b + l) 



In the limit w ^ pc equation (111 gives w oc p oc r 



and therefore ph = {w/pc )F scales as 
Ph (X , S = 6(2 — s) — s . 



(61) 



For model B2H with b 3/2 and s = 4/3 this yields ph oc 
r~^^^. Hence pfi is expected to decrease and pm = p — p.h 
to increase along the fleld lines, in agreement with what is 
observed in the simulation. 

Similar behaviour has been found in the self-similar 



solutions of Vlahakis & Konigl (2003b I, but only in cases 



where the flow is super-Alfvenic from the start (see also 
[Vlahakis at al.||2003[). In th eir trans- Alfvenic, hot-flow so- 
lutions (Vlahakis & Konigl 2003a I, pm remained constant 



throughout the thermal acceleration phase. This could be 
understood from the fact that these solutions corresponded 
to fe ~ 2 and therefore to 5 « in equation ( 61 1, resultiiig in 
constant pn and pm in the thermal acceleration regionFjln 
contrast, the super-Alfvenic solutions presented in |Vlahakis| 
& Konigl (2003b I corresponded to fe « 3/2 and hence to 



5 « — 1/3 (the same values as in our models B and D), and 
therefore they exhibited the same behaviour in the thermal 
acceleration zone as our simulated flows. 

The increase of the Foynting-to-mass flux ratio pm in 
the thermal acceleration regime leads to a rather unusual be- 
haviour of the azimuthal velocity. The right panel of Fig. [22] 
shows the variation of rv'^ along the same magnetic surface 
in models B2 and B2H. For the cold jet it always grows 



^ As was shown analytically in the magnetodynamic self-similar 
solutions of [Narayan et alT] ( |2007| l, the field-line shape is 
z oc j-2/(2— F;^ where -F is a constant parameter entering the 
self-similarity expression of the magnetic fiux function, ^ = 
T(r / z). The MHD self-similar solutions follow the same scaling 
in their force- free regime. The trans- Alfvenic solutions presented 
in | Vlahakis &: Konigl|p003a^ were characterized by _F 1, which 
implies 6 2. Note in this connection that the F = \ magneto- 
dynamic solution is exactly the paraboloidal force-free solution 



presented by Blandford 1 1976 I 
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Figure 21. The ratio of flow half-angle, 8^, to the Mach angle, Sm, across the jet for models A (left panel), B3 (middle panel) and F 
(right panel). For model A the depicted cross-sections are at i? = 10 (solid hne), R = 10^ (dashed line), R = 10^ (dash-dotted hne), 
R = 10* (dotted line) and R = 10^ (dash-triple-dotted line). For model B3 the depicted cross-sections are at r; = 5 X 10'^ (thin solid line), 
ri = 5x 10^ (dashed line), r; = 5 X 10* (dash-dotted line), = 5 X 10^ (dotted Une), r] = 5 X 10^ (dash-triple-dotted line) and r; = 5 X 10^ 
(thick solid line). For model F the depicted cross-sections are at r; = 1.5 X 10'^ (thin solid line), r) = 5 X 10'^ (dashed line), rj = 1.5 X 10* 
(dash-dotted line), and rj = 1.5 X 10^ (dotted line). The curves in the right panel dive to zero when the flow becomes sub-magnetosonic. 




with cylindrical radius and hence with the distance from 
the jet origin. This reflects the fact that the plasma is be- 
ing spun up by the rotating ma gne tic field: in this case 
\Bp/B^\ > \vp/v^\ in equation (14l and fa rQ x r. 
However, in the hot jet ru*^ (and therefore also w'*) initially 
decreases with increasing r and even attains negative val- 
ues, indicating counter rotation of the plasma. Eventually 
the cold-jet behaviour is restored, with the switch taking 
place at the turning point of /irn- The decrease in rv'^ when 
/im increases along a field line follows from the following 
relation. 



1 - in/fic^ 
1 — Mm/M 



(62) 



obtained by combining equations ( 15 I and | |l6[ )P] Physically, 
the increase in /i,„ implies that the magnetic contribution 
to the total angular momentum per unit rest mass goes up 



(see equations 20 and 15 1, which, by the conservation of I 



along a field line (and taking account of energy conservation) 
implies that the specific material angular momentum rv"^ 
must decline. 

The efficiency of the acceleration in model B2H is higher 
than in the cold models, as can be seen in Fig. [13] This is 



* The inequality IQ/fic^ < 1 always holds in trans- Alfvenic flows, 
since (Ifl/^c^)^/'^ equals the value of r/ric at the Alfven surface 
(e.g. [Vlahakis fc Komgl|2003a[ l, and the Alfven surface is located 
closer to the source than the light cylinder (with the two surfaces 
almost coinciding for highly magnetized flows). 
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connected to the behaviour of the function S. The increase of 
S during the thermal acceleration phase results in a higher 
iSf, the value of the function cS at the fast magnetosonic 
surface. In addition, the asymptotic value 5oo is smaller than 



interior field lines would tend to cylinders (see Chiueh et al 



in cold models (see Fig. 22 l. Both effects result in a higher 
value of Foo/m (see equation 31 1. 



5.6 Comparison with semi-analytic solutions 

As discussed in Section [T] it is possible to find exact so- 
lutions of the relativistic MHD equations by assuming ra- 
dial self-similarity ( Li et al.||1992 Contopoulos 1994 Vla- 
hakis & Konigl 2003a|b 20041. Due to the mathematical 
complexity of the equations, these are the only possible 
exact semi-analytic solutions describing cold or polytropic 
flows (Vlahakis & Konigl 2003a I. Similarly to their non- 



relativistic counterparts (the Blandford-Payne-type mod- 
els), they successfully capture the physics of magnetically 
driven jets and yield the general characteristics of the flow 
acceleration and collimation|^ In particular, the results of 
Vlahakis & Konigl ( 2003a|b] ) for ultra-relativistic GRB jets 
follow the general scaling relationships derived here. In fact, 



the scaling F oc 



corresponding to a streamline shape 



r (for 1 < fe < 2; the regime (iii) of equation 371, was 



first presented in Vlahakis & Konigl 1 2003b I. Note in this 



connection that both the F ~ r/ric (equation 38 1 and the 
F ~ z/r (equation |41[ ) scalings exhibited by our solutions 
could be captured through the basic radial-self-similarity 
ansatz F = T{6) because both r/ric and z/ric are functions 
of the polar angle 9 in the self-similar solutions. The semi- 
analytic solutions exhibit as high an acceleration efficiency 
( > 50%) as the simulated b < 2 solutions, and, correspond- 
ingly, have a similar value for the asymptotic shape function 
(5oo ~ 1/2; Vlahakisl|2004b I . Self-coUimation also acts in a 
similar way in both types of solution, with the inner field 
lines at any given height z being better aligned with the 
rotation axis than the poloidal field at larger values of r. 

Despite their qualitative similarity in regard to the ac- 
celeration and coUimation processes, the semi-analytic and 
numerical solutions do of course differ in their details, re- 
flecting the fact that in the self-similar model the angular 
velocity at the base necessarily scales as 1/r and that only 
one current-flow regime is allowed. In particular, the spatial 
distributions of the integrals of motion is not the same in 
these two cases. For example, the energy integral, which is 
constant in the self-similar model, is roughly proportional 
to the magnetic flux function in the simulated uniform- 
rotation jets, and the adiabat Q, which is given as a power 
of the magnetic flux function in the self-similar model, is 
a global constant in the simulations. We also note that, 
while the far-asymptotic (beyond the acceleration region) 
flow shape in the self-similar models is either cylindrical or 
conical, only the innermost field lines become cylindrical in 
the simulated jets, whereas further out the streamlines re- 
main paraboloidal. However, this is evidently related to the 
imposed boundary shape, and we can expect that, if the 
flow were followed to still larger distances, even more of the 



^ The self-similar solutions of |Vlahakis fc Konigl] | |2003a^ have a 
line-shape z oc (see footnote [7| and thus most closely resemble 
our model C. 



1991 1 or (in the case of an initially "hot" flow) to cones. 



The high acceleration efficiency inferred from the self- 
similar and numerical solutions for non-radial, relativistic 
MHD outflows was also deduced by |Beskin fc Nokhrina| 
( 2006 1 on the basis of a perturbative analysis around a 
parabolic (2 cx r^) flow. These authors found that the 
Lorentz factor increases with distance from the origin as F oc 
z^^'^, in agreement with our general result for paraboloidal 
jets of this type, F ~ z/r. 



6 APPLICATION TO GRB JETS 

The observational study of GRBs has not yet reached the 
stage where the basic parameters of the flows producing 
prompt 7-ray emission and afterglows have become well es- 
tablished. There is no general consensus yet on the angular 
structure, degree of coUimation, distance from the central 
source or composition of GRB jets. These parameters may 
vary significantly from burst to burst. The anisotropy of 7- 
ray emission due to relativistic beaming further complicates 
the problem as the same burst could have a very different 
appearance when observed from different viewing angles. In 
this section we test our theory against the current, not yet 
very stringent, observational constraints and provide a guide 
for future observations. 

The maximum terminal Lorentz factors in our numeri- 
cal models of parabolic jets, ~ 100 — 300, are close to those 
inferred for long/soft GRB jets and also high enough to en- 
sure that we have captured the properties of magnetic accel- 
eration in the ultra-relativistic regime. Although real GRB 
jets may be even faster (e.g. Lithwick fc Sari|2001 1, the an- 
alytic results verified by our numerical study can be applied 
to such jets with a high degree of confidence. 

To make detailed comparisons between our theory and 
the observations we need to determine the characteristic 
light-cylinder radius at the source of the jets. In the case 
of a millisecond magnetar 



cT 
2^ 



5 X 10' 



Vims/ 



and for a maximally rotating black hole 
M 

Thus, L — lO^cm is a suitable reference length-scale for this 
application. 

Given the extended nature of magnetic acceleration, the 
first question that one has to address is whether the Lorentz 
factors deduced from observations can be reached in our 
model on the inferred scale of the 7-ray emission region. 



According to equations ( 38 1 and ( 42 1 



for paraboloidal jets with & = 3/2 and 6 = 3, and 



for paraboloidal jets with b — 2. These estimates are lower 
than the distance to the 7-ray production region inferred 
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from the burst variability in the internal-shocks model of 
GRBs, 



r^cSt = 3 X 10^ 



St 



Vioo; vo.is 



where 5t is the internal variability time-scale (e.g. |Piran 
20051. In fact, recent Swift obser vations indicate even larger 
distances (~ lO^'^ - 10^'' cm; e.g. |Lyutikov|2006a Kumar et 
|al.|2007| ). The theory thus appears to be consistent with the 
observations in this respect. 

We emphasize that the above results have been derived 
in the context of ideal and axisymmetric MHD. In real- 
ity, various instabilities, and in particular non-axisymmetric, 
current-driven ones occurring near the jet axis, may result 
in magnetic reconnection and dissipation. It is interesting to 
note in this connection that the dissipation of Poynting flux 
would naturally generate a negative magnetic pressure gra- 
dient (associated with the azimuthal fleld component) along 
the flow and that this process was argued to be capable, on 
its own, to accelerate the flow to a high Lorentz factor (e.g. 
[Dren khahn fc Spruit|2002||Drenkhahn|2002| |. In this respect 
our ideal-MHD simulations may be yielding only lower limits 
on the terminal Lorentz factor in the modelled jets. 

A related issue is whether there is an adequate conflning 
medium, as required for the establishment of the "power- 
law" acceleration regime described by equation ( 42 I . If the 



confinement of a long/soft GRB jet is provided only by 
the envelope of the progenitor massive star, as proposed 



by Tchekhovskoy et al. (20081, the acceleration would need 



to take place on a scale smaller than the stellar radius, 
~ 10" - 10^^ cm. Downstream of the stellar surface the 
jet is expected to enter the regime of "free" (ballistic) ex- 
pansion, as in our model E, which is characterized by a less 
efficient magnetic acceleration]^ But even this rather re- 
strictive constraint on the size of the acceleration region, 
and hence on Too, is in principle consistent with the the- 
ory. An alternative possibility is that the GRB outflow is 
confined by a wind launched from the surface of a disc that 
surrounds the central object (e.g. Levinson fc Eichler|2000 1. 
This mechanism is a prime candidate for the confinement 
of short/hard GRB outflows, which evidently do not orig- 
inate inside a star. In this case the coUimation might be 
attained smoothly, with the disc-driven and central object- 
driven components constituting parts of a coherent outflow 
configuration (e.g. Tchekhovskoy et al. [2008 1. However, the 
outflow may also involve shocks formed at the interface of 
these two components (e.g. Bromberg fc Levinson|[2007 1 . If 
the GRB jet and disc outflow commence at the same time. 



It has been suggested that matter-dominated GRB jets could 
remain confined by the expanding cocoon of rclativistically hot 
shocked jet material after they break out through the stellar sur- 
face (e.g. Ram irez-Ruiz et al.|20 02 ) and could continue to acceler- 
ate during that phase (e.g. |Lazzati fc Bcgelman 2005j . In contrast, 
Poynting-dominated jets do not inflate large cocoons but instead 
create the so-called "nose cones" (e.g. |Komissarov|1999b[ l. In fact, 
given the low compression ratio of a fast shock in a magnetically 
dominated plasma, a jet termination shock is unlikely to form be- 
fore the jet emerges from the star — instead, the jet would have 
the form of a super- Alfvenic but sub— fast-magnetosonic outflow, 
as has been observed in recent computer simulations (e.g. |Komis-| 
sarov fc Barkov|2007| [Barkov fc Komissarov|2008| . 



the spatial extent of the confining medium in this picture 
can be estimated as 



7?wind ~ 3 X 10" 



9 / f wind 



O.lc 



where «wind is the mean wind speed over this distance and At 
is the GRB duration (normalized here to a fiducial value ap- 
propriate for a short/hard burst). This should be compared 
with the above theoretical relationships between R and F, 
which for F = 30 (a fid ucial value for the lower limit on Foo in 
. Nakar||2007[ ) yields _R « 3 x 10^° cm 
3 and 7? « 9 X 10** cm for b = 2. 



short/hard GRBs; e.g. 
for 6 = 3/2 or 6 = 



This comparison indicates that, over the time At, a mod- 
erately relativistic disc outfiow could form a sheath around 
the jet acceleration region. Given that the size of a disc that 
forms during a binary (NS-NS or NS-BH) merger that gives 
rise to a short/hard GRB event is not expected to exceed 
a few times 10^ cm (i.e. significantly less than than the ex- 
pected cylindrical radius of the jet in the main acceleration 
region), meaningful confinement would be attained only if 
the wind had sufficiently large inertia, which would require 



the wind-to-jet total energy ratio to be ^ 1 (cf. Levinson 
|fc Eichler|2000| ). If the initial magnetizations of short /hard 
and long/soft GRB outflows are comparable, this scenario 
provides a plausible explanation of the flnding (from the 
best available current data) that short-GRB jets are on av- 
erage less relativistic than their long-duration counterparts. 
A concomitant prediction, which could be tested when more 
afterglow data for short/hard GRBs become available, is 
that short/hard GRB outflows should also be less well col- 
limated, on average, than long/soft ones. 

The internal-shocks model envisions the prompt GRB 
emission to be powered by the collision of successively 
ejected relativistic "shells" (e.g. |Piran||2005| ). This scenario 
requires the jet to be kinetic 

energy-dominated on the scale of the emission region; 
otherwise, the flow deceleration and dissipation at fast 
shocks is too weak (or else, if the flow is inhomogeneous, the 
energy requirements are strongly increased). The numeri- 
cal solutions presented in this paper have demonstrated the 
possibility of efficient conversion of Poynting flux into bulk 
kinetic energy, with > 50% efficiency attained by the end 
of the power-law-like acceleration regime. However, the dis- 
tance Ry of the prompt emission region from the central 
source imposes a constraint on the initial magnetization of 



GRB jets in this model. Using equation (451, we obtain 



M ~ 2Foo < 



2('^7/nc)''- 

2(r^/nc) 



if 
if 



b<2 
b>2 



(63) 



For paraboloidal jets with b = 3/2 or 6 = 3 this gives (setting 

-Rlc ~ He) 

M< 430 (-^) , 
V lO^^cm/ 

whereas for & = 2 we obtain 

By approximating VacMjC? ~ £-1 At, where E the outflow 
kinetic energy as inferred from afterglow observations and 
At is the burst duration, we estimate the mass outflow rate 
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in the jet to be 

Mi ^ 5.6 X 10~' I — ^, I I I I — TT I ivip) 

' ' IQSierg J VlOs/ VlOV " 

where we normalized by values appropriate to long/soft 
bursts. This is very much lower than the expected mass 
accretion rate onto the central black hole in the coUa psar 
model (~ 0.05 - 1 Mqs'^ e.g. |Popham et al.||l999[ ) and 
constitutes the so-called "baryon loading problem" in GRB 
source models. Such a comparatively low mass outflow rate 
might be produced if the GRB-emitting outflow originates 
on magnetic field lines that thread the horizon of a spinning 
black hole and tap its rotational energy via the Blandford- 
Znajek mechanism (e.g. [Levinson fc Eichler||1993l ); in this 
case the flow would initially be baryon-free and would re- 
quire a baryon-injection mechanism as it propagates out- 
ward. Alternatively, jets launched from an accretion disc 
may experience such a low mass loading if they are ini- 
tially thermally driven along magnetic fleld lines inclined 
at a small (^15°) angle to the rotation axis (Barzilay & 

Levinson|2008| n 



portionality between the Lorentz factor along a poloidal 
magnetic surface and tan 6^ for that surface for paraboloidal 



The internal-shocks model of GRBs has been questioned 
on account of the relatively high emission efficiency that 
it requires, and these challenges have become significantly 
stronger following observations made by Swift (e.g. Gra- 
not et al. 2006 Kumar et al. 20071. Various suggestions 



have been made (and continue to be made) in the litera- 
ture for reconciling this scenario with the observations (e.g. 
Kobayashi fc Zhang|2007 1 or else for modifying or replacing 



it. Perhaps the main alternative picture proposed to date 
is based on the assumption that the prompt high-energy 
emission is produced directly from the dissipation of mag- 
netic energy without requiring it to be converted into ki- 
netic energy first (e.g. Kumar et al. 2007 1, which circumvents 
the efficiency problem that has troubled the internal-shocks 
model. Although magnetic dissipation could in principle oc- 
cur also in the context of the MHD model 



(e.g. 



Drenkhahn 



fc Spruit|2002 1, perhaps the most extreme realization of this 
idea occurs within the framework of the magnetodynamics 
scenario, in which GRB outflows are regarded as remain- 
ing Poynting flux-dominated (and sub-fast-magnetosonic) 
in the 7-ray emission region (e.g. |Blandford|2002| |Lyutikov| 
2006b I. In this scenario, neither the internal nor the reverse 



shocks of the standard model would develop, which could be 
the basis for an observational test|3 

As we discussed in Section [5. 2| a key prediction of the 
magnetic acceleration model is the approximate inverse pro- 



11 It was also proposed that the problem could be alleviated in 
a magnetically driven disc outflow that is initially neutron rich 
and hot if the neutrons decouple from the protons well before 
the latter attain their terminal Lorentz factor (see |Vlahakis at al.| 
[2003 and Fuller ct al. 2000) . There are indications from studies 
of discs around non-rotating black holes that this might not work 
in practice because outflows may be required to be comparatively 
massive to remain neutron rich (e.g. |Levinson||2006| [Barzilay &:| 
[Levinson][2008[ l, but this conclusion still needs to be verified in 
the case of discs around rapidly rotating black holes. 
1^ Note in this connection that, in some of the proposed interpre- 
tations of the Swift data (e.g. [Uhm fc Beloborodov|2007| [Genet et| 
[al. |2007| , the entire afterglow emission is attributed to a reverse 
shock that is driven into the ejecta. 



jets with 1 < 6 < 2 (see equation 43 1. For a small opening an- 
gle and b not very close to 1 this result can be approximated 
as F^v ~ 1. This implies that GRB outflows with b < 2 that 
attain F ~ 100, the approximate inferred lower limit for 
long/soft GRBs, must have 6v ~ 0.6°, essentially indepen- 
dent of the details of the acceleration process. When b > 2, 
rOv ~ b~^{R/ric)~^^~^^''^ decreases with R in the magnetic 
acceleration region, implying an even smaller value of 9^ at 
the end of this zone. The relation F^v ~ 1 may be useful 
for differentiating between magnetic and fireball models of 
GRB flows. Indeed, this property is generic to the magnetic 
acceleration mechanism, whereas for the thermal acceler- 
ation the terminal bulk Lorentz factor is essentially given 
by the thermal Lorentz factor at the base of the flow and 
is fairly independent on the flow coUimation, which means 
that the product FO^ can in principle become ^ 1. Interest- 
ingly, one of the proposals made for interpreting the appar- 
ent GRB "tails" observed by Swift invokes a GRB-emitting 
outflow component whose opening half-angle must be < 1° 
(Panaitcscu 2007) . While the currently available data are 
not sufficient for favouring this interpretation over other 
suggested explanations of the "tails," it is noteworthy that 



the requirement arrived at by Panaitescu ( 2007 1 on strictly 



phenomenological grounds is consistent with a distinguish- 
ing property of the magnetic acceleration model. It is also 
noteworthy that there is already at least one source (GRB 
070401) in which such a small opening half-angle has been 
inferred directly from a measurement of an early break in the 
X-ray afterglow light curve ( Kamble et al.|2008 1. Such small 
asymptotic opening angles and even F^v ^ 1 could in prin- 
ciple be attained also in purely hydrodynamical jet models, 
although this would require a very high efficiency of coUima- 
tion and acceleration within the stellar interior. Specifically, 
the jets would need to emerge from the star with 6^ < 1° 
and r>l/9v — 60, which, in view of recent analytic and 
numerical studies (e.g. Lazzati fc Begelman||2005 Morsony 
[et al. |2007[ ), is unlikely to be achieved in practice. 

The original fireball model for GRB jets envisions a uni- 
form conical outflow that becomes accelerated to Lorentz 
factors F 3> l/6v and predicts that during the afterglow 
phase the Lorentz factor of the forward shock driven by 
the jet into the ambient medium will decrease to values 
< l/8v The observational consequence of this transition is 
a panchromatic break in the afterglow light curve (referred 
to as the "jet break") occurring when 6vF becomes ~ 1 (e.g. 
Rhoads 1999 Sari et al. 19991. In view of the results pre- 



sented in this paper, the predictions of the MHD model for 
GRB outflows that are efficiently accelerated — and there- 
fore necessarily confined (by either thermal, magnetic or ram 
pressure) during the acceleration phase — are radically dif- 
ferent. Specifically, the MHD model predicts that the after- 
glow light curve would exhibit either a very early jet break 
(in cases where F^v ~ 1 at the end of the acceleration phase, 
as expected in jets with 6 < 2) or no jet break at all (if 
F^v < 1 at the end of the magnetic acceleration region, as 
expected in jets with b > 2) This prediction is seemingly 



1'^ If the low current detection rate of jet breaks in the early 
afterglow light curves of GRB sources would prove to be more 
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at odds with the inference from a number of pre-Swift GRB 
sources of breaks of this type occurring on a time-scale of 
days (see e.g. Liang & Zhang 2005 for a compilation). The 
paucity of "textbook" jet breaks in Swift GRB sources (e.g. 
Liang et al.|2008 l, which has even cast doubts on the inter- 



the kinetic power across the jet directly affects the evolu- 



tion of the light curve of the GRB afterglow (e.g. Granot 



2005 1 as well as the statistical properties of a GRB sample 



pretation of the alleged pie-Swift jet breaks, points to one 
way out of this dilemma: it may be that indeed there are no 
bona fide jet breaks at later times. We recall, however, that 
the jet-break interpretation lies at the basis of the identi- 
fication of GRB outflows as coUimated jets, which has sig- 
nificantly reduced the otherwise prohibitive energy require- 
ments in some sources. Alternatively, it could be that the 
difficulties in finding late-time jet breaks in Swift sources are 
to a large extent observational (e.g. |Zhang„2007| ), in which 
case other explanations for late-break candidates must be 
sought. 

One natural possibility is that the outfiow possesses 
more than one kinematic component. In its simplest incarna- 
tion, this is the "two component" model, which envisions the 
prompt emission to originate in an ultra-relativistic, highly 
collimated jet and the afterglow emission to be dominated 
by a less relativistic, wider outflow component. The sugges- 
tion in Panaitescu ( 2007 1 and in Kamble et al. ( 2008 1 that 



the 7-ray emitting jet is very narrow was made in the con- 
text of this model, and a similar picture was used by lGranot] 



et al. ( 2006 1 to explain other aspects of the early GRB X-ray 
emission measured by Swift (see also Zhang 2007). In fact, a 



two-component outflow conflguration had already been pro- 
posed in the pie-Swift era to account for certain observations 
(e.g. Berger et al.||2003b l and as a means of alleviating the 
efficiency requirements on the internal-shocks model ( |Peng| 
|et al.,2005j . The separation into two components could arise 
either from an interaction of the outfiow with the envelope of 
a massive progenitor star or represent an intrinsic property 
of the central engine (see |Peng et al.|2005| for a summary of 
some specific proposals). In the context of the magnetically 
driven outfiow model, there are at least two possibilities for 
an intrinsic origin. First, neutron-rich, hot outfiow may split 
into two components when the neutrons and protons decou- 
ple before the protons have attained their terminal Lorentz 
factor (Vlahakis at al.|;2003). Second, a baryon-poor ultra- 



relativistic outfiow launched from the black hole can be sur- 
rounded by a magnetically driven, relativistic outfiow from 
the accretion disc itself (see Granot et al.|2006|p^ We stress 
that, in reality, the outfiow may be more complex than in 
the schematic "two component" picture sketched above. For 
example, inhomogeneities in the accretion fiow may result in 
several distinct outfiow components emerging from the disc, 
associated, perhaps, with isolated magnetic fiux tubes that 
thread the disc at different locations. Phenomenologically, 
this situation might resemble the "patchy shell" scenario 
considered by Kumar & Piran ( 2000 1 . 



The distribution of the terminal Lorentz factor and of 



than just the result of observational difficulties, this could be 
an indication, when interpreted in the context of the magnetic 
acceleration model, that these jets are characterized by effective 
shape-function exponents b > 2. 

In the latter scenario, the disc wind could provide a ready 
source for seeding the central funnel with baryons (e.g. |Levinson| 
I&: Eichlcr 2003| and could also help coUimate the interior outflow 
( [Levinson fc Eichler|2000 i. 



(e.g. Nakar et al. 2004 1 and the detectability of "orphan" 
afterglows (afterglows detected without an associated GRB; 
e.g. Nakar Sz Piran|2003 1. One could in turn attempt to use 
such observations to probe the jet structure and to test the 
underlying acceleration and coUimation models. With this 
in mind, we present in Fig. [23] illustrative asymptotic distri- 
butions of the Lorentz factor and of the kinetic power from 
our simulations. 

We consider first the Too distribution. The top left panel 
of Fig. [23] shows that in all of the cases the Lorentz factor 
decreases toward the axis — this is a generic feature of the 
axisymmetric, ideal-MHD acceleration mechanism as the az- 
imuthal magnetic field and hence the Poynting flux vanish 
along the symmetry axis. This feature may not, however, 
be as pronounced when non-axisymmetric instabilities and 
resistive dissipation of magnetic energy (which are not in- 
corporated into our study) are taken into account. In fact, 
we find that even in our solutions F 7^ 1 at S = because of 
numerical dissipation. In the case of an initially hot outfiow 
r{9 = 0) > 1 is due to the thermal acceleration. In initially 
cold outfiows that have uniform rotation and mass density 
distribution at the base F peaks at the jet boundary. It is 
seen, however, that if the fiow is initially hot the anisotropy 
of the Lorentz factor distribution within the jet is reduced. 
Uniform rotation is a robust prediction of models with a 
magnetar or a magnetized black hole as a central rotator. 
The assumption a uniform mass-fiux distribution at the jet 
base is more of an approximation: for example, when the 
central source is a black hole the degree of baryon loading 
is likely to be higher near the jet boundary due to vari- 
ous boundary interactions with the jet surroundings (e.g. 
,Lcvinson fc Eichlcr 2003). Such a mass distribution would 
lead to lower terminal Lorentz factors near the boundary 
compared to that found in our simulations. If the inner re- 
gions of an accretion disc contribute to the magnetic driving 
of the GRB-emitting outfiow component then a model with 
differential rotation, with Q, decreasing away from the cen- 
tre, is more suitable. As seen from the figure, in this case 
the terminal Lorentz factor peaks at intermediate angles. 
In practice it may, however, be difficult to distinguish this 
case from that of uniform rotation with nonuniform mass 
loading. 

Turning now to the distribution of energy fiux across 
the jets in the asymptotic regime, we recall that the obser- 
vational consequences of this energy are strongly influenced 
by relativistic beaming — whenever a fraction of this energy 
is dissipated and converted into radiation, this radiation will 
be beamed in the direction of motion of the corresponding 
fluid element, given by 9^. Most phenomenological models 
of GRBs have assumed that the jet is conical and has ra- 
dial streamlines. Thus, the streamline angle, 9^, is equal to 
9, the polar angle of the fluid element. In our model the 
streamlines are curved and asymptotically their shape is 
close to that of the boundary (with the exception of the 
cylindrical core). Hence we have 9^ — 9/a. Consider a sur- 
face element normal to the 77 coordinate lines (streamlines) 
dSr, = ^/g4^^pg^d(j}d^, where gr00 and (755 are components of 
the metric tensor. Since in the asymptotic regime 9, 9^ <C I, 
we can write 



"^/"d^v, g^,/, = and = z^^" (see 
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Figure 23. Angular distributions of the Lorentz factor (top left panel), the kinetic power per unit solid angle in the local direction of 
the flow (e, top right panel), the kinetic power per annulus of unit angular size bottom left panel) and ed^ (bottom right panel) in 
the asymptotic regime, plotted as functions of polar angle. The variable e is given in units of cB^L^ /An, and when it is multiplied by 8 
or 6^ the polar angle is measured in radians. Note, however, that the polar angle along the horizontal axis is given in degrees. The solid 
lines show model B2, the dashed lines model B2H, and the dash-dotted lines model D. 



Appendix A of Paper I) , and hence 



where dui — O^dOvdcp is the soUd angle defined by the tan- 
gents to the streamlines passing through the surface element. 
The power per unit solid angle is then given by 

dC/did = e — S^a^R^ , 

where 5"' is the component of the energy flux density in the 
rj direction. 

The top right panel of Fig. |23] shows the distribution 
of kinetic power per unit solid angle, e (the total power has 
a very similar distribution). One can see that in all models 
it peaks at, or very close to, ^ = 0. The reason for this 
behaviour, which seemingly conflicts with the Lorentz factor 
distribution shown in the left panel of the figure, is that the 
density distribution across the jet is highly nonuniform, with 
the mass density strongly peaking near the symmetry axis 
on account of the enhanced coUimation of the fiow in that 
region (see Figs.[2]|5|. The bottom left panel of Fig. [23| shows 
the distribution of eO: this quantity tells us how the jet power 
is distributed between annuli of equal size in 6. One can see 
that in model D (differential rotation) more power comes 
from the intermediate annuli, in model B (uniform rotation 



at the base) from the outer annuli, and that a significant core 
component emerges in model B2H (initially hot jet). Note, 
however, that the distributions of the Lorentz factor and the 
power depend on the choices of the density, magnetic fiux 
and angular velocity distributions at the inlet boundary, so 
different profiles may be possible. 

The derived distributions of r(^) and e[9) are markedly 
different from those commonly adopted in phenomenologi- 
cal GRB jet models, which either take them to be uniform 
within the jet half-opening angle 9j or else assume that the 
ffow has a universal structure, with e bei ng a Gaussian or 



Rossi et al. 



2002 



a power-law in 6 (in particular 
— compare with the bottom right panel of Fig. 23 1 out- 
side a uniform-core region]^ The structure exhibited by our 



In the force-free electromagnetic model for GRBs it is en- 
visioned that the current flows along the axis of rotation and 
returns through the equatorial plane; this yields an energ y distri- 
bution oc in the associated electromagnetic shell (e.g. |Bland-| 
|ford|2002] |Lyutikov|2006b"l l . A universal structured outflow with 
e oc could potentially also be produced when a relativis- 

tic GRB jet with possibly a different initial energy distribution 
breaks out through the surface of a massive progenitor star | |Laz-| 
zati fc Begelman|20"05 l. 
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model jets is also different from that of a "hollow cone," 
where the flow occupies the region £ [9j — AO, 6j] (e.g. 



Eichler & Levinson 2004 Lazzati & Begelman 2005). Al- 
though the distribution of Lorentz factors is reminiscent of 
such a cone, the distribution of kinetic power actually peaks 
near the symmetry axis. Moreover, in contrast with the phe- 
nomenological hollow-cone models considered in the litera- 
ture, in which AO <^ 9j, our solutions yield configurations 



with AO 



The detailed observational implications of 



these structures remain to be explored. 



7 CONCLUSION 

In this paper we extend our previous numerical study of 
magnetically accelerated relativistic jets (Paper I) from the 
case of terminal Lorentz factors Too ~ 10, appropriate to 
AGN jets, to Too ^ 10^, appropriate to GRB jets. The larger 
values of Too reached in the present study enable us to com- 
pare results of our simulations, carried out using the equa- 
tions of special-relativistic ideal MHD, with the asymptotic 
analytic formulae that we obtain from the constituent equa- 
tions in the limit F 3> 1. Our analysis of the results also 
benefits from a comparison with semi-analytic solutions that 
were derived under the assumption of radial self-similarity. 
We can summarize our conclusions regarding the magnetic 
acceleration of ultra-relativistic outflows as follows. 

(i) Our simulations verify that the MHD acceleration 
mechanism remains robust even when the terminal Lorentz 
factors reach the ultra-relativistic regime (roo>102). The 
simulated flows rapidly settle into quasi-steady and seem- 
ingly stable configurations. A complete model would need 
to incorporate non-axisymmetric effects, which we have not 
considered. 

(ii) A key property of magnetically driven relativistic 
flows in the ideal-MHD regime is the spatially extended 
nature of their acceleration. This property, which was first 
revealed by the self-similar solutions and subsequently con- 
firmed in the moderately relativistic regime by the simula- 
tions reported in Paper I, is also a distinguishing charac- 
teristic of jets accelerated to ultra-relativistic speeds. For 
initially Poynting flux-dominated jets whose magnetic flux 
surfaces can be approximated by paraboloids of the form 
z oc r'' (with 6 > 1), the Lorentz factor during the main 
magnetic acceleration phase increases as F ~ {b/\/b — l)z/r 
when 1 < 6 < 2 and as F ~ r /r\c when fe = 1 or 6 > 2. After 
the (increasing) kinetic energy flux becomes comparable to 
the (decreasing) Poynting flux the growth of F saturates, 
and thereafter it increases at a much slower rate. (We have 
not been able to reach this phase in models with b > 2 due 
to the limitations of our numerical method.) 

(iii) The conversion efficiency Foo//i of total injected en- 
ergy to kinetic energy at the end of the power-law acceler- 
ation phase lies in the range 55 — 75% for the initially cold 
simulated paraboloidal flows whose effective exponents lie in 
the range 1 < & < 2; the efficiency is smaller the larger the 
initial magnetization (or, equivalently, the higher the value 
of Foo). A higher efficiency is attained in jets with fe < 2 
that are initially relativistically hot than in the correspond- 
ing initially-cold outflows: in this case a measurable fraction 
(> 50% in the example that we show) of the thermal energy 



flux is at first converted into Poynting fiux, thereby reduc- 
ing the initial thermal acceleration of the flow and enhancing 
the subsequent magnetic acceleration. 

(iv) In our simulations the flow is confined by a rigid wall 
whose shape is described by z oc r", with a ranging from 2/3 
to 3. We have conducted a detailed analytic investigation of 
the relationship between a confining pressure distribution 
of the form pcxt oc z~°' and the shape of the jet boundary 
in the asymptotic regime of the magnetic acceleration zone. 
We found that there is a one-to-one correspondence between 
the functional forms of the pressure distribution and of the 
boundary shape. Except for one special case (for which a re- 
mains close to 2), the jet becomes an exact paraboloid of the 
form given above, with a = 4/a > 2 for a < 2 and 1 < a < 2 
for a — 2. When a > 2 the jet cannot maintain pressure 
equilibrium with the ambient medium and asymptotes to a 
conical shape. This situation is reproduced in our simula- 
tions by unconfined flows as well as by flows with a < 1. 
In this case the outer regions of the jet become causally 
disconnected (the local opening half-angle of the field lines 
becomes larger than the local half-angle of the Mach cone 
of fast-magnetosonic waves) , and only the innermost regions 
continue to coUimate and accelerate. 

(v) We find that for all current-carrying jets (irrespective 
of whether the return current flows inside or outside the 
jet) the innermost field lines are more strongly collimated 
than the exterior ones, indicating "self collimation" by the 
magnetic hoop stress (see also Paper I). This redistribution 
of the poloidal field lines within the jet is directly responsible 
for the high acceleration efficiency of the flow. 

We have applied our results to GRB sources, taking into 
account the constraints imposed by the detected prompt and 
afterglow emission on the properties of the ultrarelativistic 
jets that evidently give rise to the GRB phenomenon. Our 
main conclusions are: 

(i) Initially Poynting flux-dominated outflows can be 
magnetically accelerated to a Lorentz factor exceeding the 
minimum (F ~ 10^) inferred in long/soft GRBs within a 
distance of ~ lO" - 10^^ cm from a rapidly rotating stellar- 
mass black hole or a millisecond magnetar. Thus, most of 
the acceleration of long/soft GRB jets can be achieved in- 
side a typical progenitor star in the coUapsar model, whose 
envelope provides a natural confining environment for the 
jets. Lack of confinement outside of the star may result in 
a radial outflow characterized by loss of causal connectiv- 
ity across the jet and inefficient acceleration. An alternative 
confinement mechanism that is of particular relevance to 
short/hard GRBs, which likely form through a merger of 
compact stars rather than in the collapse of a massive star, 
is a disc wind. The MHD acceleration mechanism implies 
that the minimum bulk Lorentz factor inferred in short /hard 
GRBs (F ~ 30) could be attained within the distance that 
such a wind covers over the burst duration if the disc out- 
flow (which might also be driven magnetically) has at least 
a moderately relativistic speed (~ 0.1 — Ic). 

(ii) The MHD acceleration model entails a high ( > 50%) 
asymptotic conversion efficiency of injected magnetic and 
thermal energy into bulk kinetic energy for effectively con- 
fined flows. If the initial magnetization is of the same or- 
der as that of the inferred Lorentz factor of a GRB jet. 



10^ 



10^ the energy conversion can be attained on 
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a spatial scale that is smaller than the indicated size of 
the prompt emission region. The model is then compati- 
ble with the internal-shocks scenario for GRBs. For a much 
higher initial magnetization the jet remains Poynting flux- 
dominated on those scales and the prompt emission has to 
be attributed to direct magnetic energy dissipation, as in the 
magnetodynamics scenario. A full treatment of the dynam- 
ics of such jets in the context of MHD would require taking 
account of the acceleration induced by the field-dissipation 
process and the use of a non-ideal, rclativistic-MIfD code. 

(iii) We have found that the MHD jet model places a 
strong constraint on the product of the Lorentz factor and 
the half-opening angle of the streamline in the asymptotic 
regime of the main acceleration region: F^v — 1 along 
paraboloidal streamlines oc r* when b < 2 (but b not too 
close to 1), and F^v oc ^"('^"^/i)) (^a,nd thus attaining even 
smaller values at the end of the main acceleration phase) 
when 6 > 2. This feature is unique to the ideal MHD mech- 
anism and could potentially serve to distinguish it from al- 
ternative models, notably the classical fireball scenario (in 
which F0V is envisioned to be » 1 at the end of the ac- 
celeration region). In particular, this property implies that, 
if long/soft GRB jets with F > 100 are magnetically accel- 
erated, they must be coUimated to 0v^l°. This result is 
consistent with one of the interpretations of the prompt 
emission "tails" discovered by Swift, although this is not 
the only possible explanation of a very small collimation an- 
gle. This relationship also indicates that the 7-ray omitting 
outflow component might exhibit a panchromatic jet break 
(corresponding to F^v — 1 decreasing from a value > 1 to a 
valure < 1) soon after it enters the afterglow phase, although 
in principle no such break need to occur (corresponding to 
cases where this product is < 1 at the end of the accelera- 
tion zone). A later jet break could potentially be seen if the 
outflow has more than one kinematic component. 

(iv) The magnetic acceleration model also makes specific 
predictions about the angular distributions of the terminal 
Lorentz factor and of the kinetic and total energy per unit 
solid angle across the jet, which can be probed by a variety 
of observations. These distributions depend on the magne- 
tization profile and the thermal energy content of the jet at 
the inlet boundary, which could in principle be constrained 
by the observations. A general characteristic of this model 
is that FaoiO) decreases with decreasing polar angle G near 
the symmetry axis. 

Although our analytic scalings have been derived in the 
limit where the jet is in the force-free regime, we empha- 
size that key parameters of interest for astrophysical ap- 
plications — including the jet velocity and the magnetic- 
to-kinetic energy conversion efficiency — could have only 
been obtained within the magneto/ij/drodynamics formalism 
that we adopted and not in the magnetodynamics (or force- 
free electrodynamics) approximation adopted in other recent 
semi-analytic and numerical investigations. Another point 
worth emphasizing is that the acceleration mechanism in- 
vestigated in this paper is identical to that considered in 
paper I. Our results are consistent with the view that the 
main difference between "superluminal" AGN jets and GRB 
jets is that the latter outflows have a higher initial magne- 
tization (and possibly also a higher initial enthalpy), which 
leads to their correspondingly higher terminal Lorentz fac- 



tors. If this picture is correct, one could use observations of 
AGN and GRB sources to deduce complementary aspects 
of the same basic phenomenon. For example, one could take 
advantage of the fact that the acceleration region in AGN 
jets is potentially resolvable by radio intorforometry to probe 
the details of the acceleration process; one could then con- 
sider the implications to GRB jets, which are not directly 
accessible to such observations. 
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APPENDIX A: SOLUTIONS OF 



EQUATION (521 



In Section |5.3| we considered the dependence of the jet 
boundary shape on the external pressure distribution, and 
we derived a second-order ordinary differential equation 
(equation |52[ ) that expresses this dependence in the asymp- 
totic regime of the main magnetic acceleration region for the 
case where the external pressure scales as pcxt oc z"" . For 
convenience, we reproduce this equation here, keeping the 
original notation: 



(Al) 



where C is a constant of the order of 1 (equation|49[). In Sec- 
tion|5.3|we obtained solutions for this equation after making 



a power-law ansatz for x{Z) (equation 53 1. In this appendix 
we consider general solutions of this equation without as- 
suming from the start that they have a power-law form. 



Al a>2 

One can identify three different regimes in this case. 
Al.l 2<Q<4 

When Q < 4 the ratio of the second (electromagnetic) to 
the third (centrifugal) terms on the left-hand side of equa- 
tion ( Al I diverges as Z ^ oo, and one can therefore neglect 



the centrifugal term in the asymptotic regime (as we also 



inferred in Section 5.3.1 1. Changing variables to 



y 



}{y) 



this equation can then be written as 

2' 



y 



dy^ dy 



y 



2- a 



/ = 0, 



whose solution is 

/(y) = CiJi/|2_„|(y) + C2yi/|2_<,|(y). 



(A2) 



(A3) 



(A4) 



For Q > 2 we have 1 — a/2 < 0, so the limit Z ^ oo 
corresponds to y ^ 0, in which case Jv{y) « y" and 
Yiy{y) ~ l/y"- Thus the term involving the Neumann 
function Yv{y) dominates, implying an asymptotic solution 
/(y) ~ C2/y^^''"~'^\ or, using the definition of y (equa- 
tion 



A2i, 



(A5) 



Thus, the solution is essentially paraboloidal (concave) with 
conical asymptotes. 
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A1.2 a- 



Changing the variable x{Z) to g{Z) = x/Z, equation (All 
becomes 



^ dZ^^^^dZ^ Z~^ 



(A6) 



This equation has an exact solution, g = const = C 
representing a flow that is conical from the start, Z = C^^'^x 



(as we already found in Section 5.3.1 1, and in which the elec- 
tromagnetic and centrifugal forces have comparable contri- 
butions. 



A1.3 a>4 

When a > 4 the ratio of the third (centrifugal) to the sec- 
ond (electromagnetic) terms on the left-hand side of equa- 
tion ( |A1[ ) diverges as ^ ^ oo, and one can therefore ne- 
glect the electromagnetic term in the asymptotic regime (as 
we also inferred in Section 5.3.11. Without this term, equa- 
tion (|ATJ becomes 

Multiplying by 2dx/dZ, this equation can be integrated to 
give 
-da; ^2 



/ dx\' 
\dZ) 



+ 



1 



D ' 



(A8) 



where _D is a constant of integration. Equation ( A8 1 can be 
further integrated to yield 



{Z-Zof 



D 



D2 



= 1, 



(A9) 



where Zq is another constant of integration. Equation ( A9 \ 



explicitly shows that the jet assumes a hyperboloidal shape 
in this case, with the asymptotes again being conical. 



A2 a = 2 



Changing variables to / = x/Z^^^ (as in equation A2| and 



q = \nZ, equation ( Al I becomes in this case 



(AlO) 



d^-y4 ^)^^ 'p 

One obvious solution is / =const, with f~^^^ = C — 1/4, 
oi Z ^ — 1/4 x^ . This solution is real only for C > 1/4. 

For df/dq 7^ 0, equation (AlO I can be multiplied by 
2df /dq and rewritten as 



dq 



= 0, 



(All) 



(A12) 



where i? is a constant of integration. This can be integrated 
to give 



or, with u = f , 

= (1 -4C)it^ +4£;it-4. 



±2 dq- 



du 



In 



E/2 



Eu 



C 



,V~Eu 



c < 



E 



4 



■ arctan 



EI2 

T 
4 



(A13) 



C 



2^Eu- 
-u + 



c > 



c 



The quantity inside the square root in the integrand must 
be positive in the asymptotic regime. This implies that u 
cannot tend to zero as Z — ^ 00, and so it either tends to 
a constant or to cxd. The first option is unacceptable as it 
gives q = const <^ Z = const. Thus, we must have u ^ 00 
as ^ — > 00. Since u = — x'^/Z, this means that the 
shape is Z (X x*" with 6 < 2. The dominant term inside the 
square root in the integrand when Z 00 is the first one, 
which implies C < 1/4 (and hence that the C > 1/4 solution 



branch of equation A13 is not physical). 

For C — 1/4, equation (A13l gives 



Z^I'^^XjE + E [In (Z/Zo)]^, or, keeping only the dominant 
terms, x — Z^^^{Ci + C2lnZ). In view of the requirement 



u — X jZ 



00 derived above, one must keep 



the logarithmic term in this solution (i.e. C2 7^ 0). We also 
keep the other term in this solution, which may be needed 
to match to the conditions at the base of the flow (see Sec- 
tion |5.3.2| . Note that this solution could not be derived from 
the ansatz employed in Section [5. 3| since it does not have a 
pure power-law form. 

For C < 1/4, equation ( |A13 1 gives (keeping only the 
dominant terms) x = CiZ^^T^Ti-^c gjjjgg^ g^g found 
above, the exponent h in Z oc x^ must be < 2 the only 
acceptable solution is a; = C1Z3 + 2 
Summarizing, 

for C > 1/4, Z = ^C- 1/4 x^, 

for C = 1/4, X = Z^^'^iCi + C2 In Z) with C2 / 0, and 
for C < 1/4, X = Ci^s + ivT^. 

The first case (C > 1/4) represents a balance between the 
poloidal curvature, electromagnetic and centrifugal terms in 



equation ( Al I, whereas in the last two cases (C < 1/4) only 



the poloidal curvature and electromagnetic terms play a role. 



A3 a <2 



As discussed in Section [5.3.3[ in this case the poloidal cur- 
vature term — the first term on the left-hand side of equa- 
tion ( Al I — can be neglected. The power-law form for x{Z) 



given in the main text is then an exact solution of this equa- 
tion. 



